Philippa Kaur

1.Seddon, P. J., Armstrong, D. P. & Maloney, R. F. Developing the science of reintroduction biology. Conserv. Biol. 21, 303–312 (2007).PubMed 
Article 
PubMed Central 

Google Scholar 
2.Seddon, P. J., Griffiths, C., Soorae, P. & Armstrong, D. P. Reversing defaunation: Restoring species in a changing world. Science 345, 406–412 (2014).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
3.Pérez, I. et al. What is wrong with current translocations? A review and a decision-making proposal. Front. Ecol. Environ. 10, 494–501 (2012).Article 

Google Scholar 
4.Brichieri-Colombi, T. A. & Moehrenschlager, A. Alignment of threat, effort, and perceived success in North American conservation translocations. Conserv. Biol. 30, 1159–1172 (2016).PubMed 
Article 
PubMed Central 

Google Scholar 
5.Swan, K. D., Lloyd, N. A. & Moehrenschlager, A. Projecting further increases in conservation translocations: A Canadian case study. Biol. Conserv. 228, 175–182 (2018).Article 

Google Scholar 
6.Wolf, C., Griffith, B., Reed, C. & Temple, S. Avian and mammalian translocations: Update and reanalysis of 1987 survey data. Conserv. Biol 10, 1142–1154 (1996).Article 

Google Scholar 
7.Breitenmoser, U., Breitenmoser-Wursten, C., Carbyn, L. N. & Funk, S. M. Assessment of carnivore reintroductions. In Carnivore Conservation (eds Gittleman, J. L. et al.) 241–280 (Cambridge University Press and The Zoological Society of London, 2001).
Google Scholar 
8.Vandel, J. M., Sthal, P., Herrenschmidt, V. & Marboutin, E. Reintroduction of the lynx into the Vosges mountain massif: From animal survival and movements to population development. Biol. Conserv. 131, 370–385 (2006).Article 

Google Scholar 
9.Saltz, D., Rowen, M. & Rubenstein, D. I. The effect of space-use patterns of reintroduced Asian wild ass on effective population size. Conserv. Biol. 14, 1852–1861 (2000).
Google Scholar 
10.Reynolds, M. H. et al. Space use and reintroduction of Laysan teal. Anim. Conserv. 15, 305–317 (2012).Article 

Google Scholar 
11.Gilpin, M. E. & Soulé, M. E. Minimum viable populations: processes of species extinction. In Conservation biology: The science of scarcity and diversity (ed. Soulé, M. E.) 19–34 (Sinauer Associates, 1986).
Google Scholar 
12.Snyder, N. F. R. et al. Limitations of captive breeding in endangered species recovery. Conserv. Biol 10, 338–348 (1996).Article 

Google Scholar 
13.McCarthy, M. A., Armstrong, D. P. & Runge, M. C. Adaptive management of reintroduction. In Reintroduction Biology (eds Ewen, J. G. et al.) 256–289 (Wiley, 2012).Chapter 

Google Scholar 
14.Colomer, M. A., Oliva-Vidal, P., Jiménez, J., Martínez, J. M. & Margalida, A. Prioritizing removal actions for the reintroduction of endangered species: Insights from bearded vulture simulation modeling. Anim. Conserv. 23, 396–406 (2020).Article 

Google Scholar 
15.Serrouya, R. et al. Saving endangered species using adaptive management. Proc. Natl. Acad. Sci. 116, 6181–6186 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
16.Rodríguez, A. & Delibes, M. El lince ibérico (Lynx pardina) en España: Distribución y problemas de conservación (Colección técnica. ICONA, 1990).
Google Scholar 
17.Gil-Sánchez, J. M. & McCain, E. B. Former range and decline of the Iberian lynx (Lynx pardinus) reconstructed using verified records. J. Mamm. 92, 1081–1090 (2011).Article 

Google Scholar 
18.Casas-Marce, M. et al. Spatiotemporal dynamics of genetic variation in the Iberian lynx long its path to extinction reconstructed with ancient DNA. Mol. Biol. Evol. 34, 2893–2907 (2017).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
19.Guzmán, J. N. El lince ibérico (Lynx pardinus) en España y Portugal (Censo-diagnóstico de sus poblaciones. DGCN, Ministerio de Medio Ambiente, 2004).
Google Scholar 
20.Simón, M. A. et al. Reverse of the decline of the endangered Iberian Lynx: Saving the Iberian Lynx. Conserv. Biol. 26, 731–736 (2012).PubMed 
Article 
PubMed Central 

Google Scholar 
21.Vargas, A. et al. The Iberian Lynx Lynx pardinus conservation breedingprogram: Iberian lynx conservation breeding program. Int. Zoo Yearbook 42, 190–198 (2008).Article 

Google Scholar 
22.Simón, M. A. et al. Recuperación del lince ibérico en España y Portugal: situación actual de sus poblaciones (XIII Congreso SECEM. Guadalajara, 2017).
Google Scholar 
23.Aranda, et al. La reintroducción del lince ibérico en Castilla-La Mancha (XIII Congreso SECEM, 2017).
Google Scholar 
24.Buderman, F., Hooten, M. B., Ivan, J. S. & Shenk, T. M. Large-scale movement behavior in a reintroduced predator population. Ecography 41, 126–139 (2017).Article 

Google Scholar 
25.Revilla, E., Wiegand, T., Palomares, F., Ferreras, P. & Delibes, M. Effects of matrix heterogeneity on animal dispersal: From individual behavior to metapopulation-level parameters. Am. Nat. 164, E130–E153 (2004).PubMed 
Article 
PubMed Central 

Google Scholar 
26.Ferreras, P. et al. Proximate and ultimate causes of dispersal in the Iberian Lynx Lynx pardinus. Behav. Ecol. 15, 31–40 (2004).Article 

Google Scholar 
27.Palomares, F. et al. Iberian Lynx in a fragmented landscape: Predispersal, dispersal, and postdispersal habitats. Conserv. Biol. 14, 809–818 (2000).Article 

Google Scholar 
28.Gastón, A. et al. Response to agriculture by a woodland species depends on cover type and behavioural state: Insights from resident and dispersing Iberian Lynx. J. Appl. Ecol. 53, 814–824 (2016).Article 

Google Scholar 
29.Blázquez-Cabrera, S. et al. Influence of separating home range and dispersal movements on characterizing corridors and effective distances. Landsc. Ecol. 31, 2355–2366 (2016).Article 

Google Scholar 
30.Ferreras, P. Landscape structure and asymmetrical inter-patch connectivity in a metapopulation of the endangered Iberian lynx. Biol. Conserv. 100, 125–136 (2001).Article 

Google Scholar 
31.Stoinski, T., Beck, B., Bloomsmith, M. & Maple, T. A. Behavioral comparison of captive-born, reintroduced golden lion tamarins and their wild-born offspring. Behaviour 140, 137–160 (2003).Article 

Google Scholar 
32.Margalida, A. et al. Uneven large-scale movement patterns in wild and reintroduced pre-adult bearded vultures: Conservation implications. PLoS One 8, e65857 (2013).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
33.Berger-TAL, O. & Saltz, D. Using the movement patterns of reintroduced animals to improve reintroduction success. Curr. Zool. 60, 515–526 (2014).Article 

Google Scholar 
34.Ferreras, P., Beltrán, J. F., Aldama, J. J. & Delibes, M. Spatial organization and land tenure system of the endangered Iberian lynx (Lynx pardinus). J. Zool. 243, 163–189 (1997).Article 

Google Scholar 
35.Benson, J. F., Chamberlain, M. J. & Leopold, B. D. Land tenure and occupation of vacant home ranges by bobcats (Lynx rufus). J. Mamm. 85, 983–988 (2004).Article 

Google Scholar 
36.López-Parra, M. et al. Change in demographic patterns of the Doñana Iberian lynx Lynx pardinus: Management implications and conservation perspectives. Oryx 46, 403–413 (2012).Article 

Google Scholar 
37.Palomares, F., Delibes, M., Revilla, E., Calzada, J. & Fedriani, J. M. Spatial ecology of Iberian Lynx and abundance of European rabbits in Southwestern Spain. Wildl. Monogr. 148, 1–136 (2001).
Google Scholar 
38.Simón, M. A. et al. Ten years conserving the Iberian lynx (Consejería de Agricultura, Pesca y Medio Ambiente. Junta de Andalucía, 2012).
Google Scholar 
39.Gautestad, A. O. Memory matters: Influence from a cognitive map on animal space use. J. Theor. Biol. 287, 26–36 (2011).PubMed 
MATH 
Article 
PubMed Central 

Google Scholar 
40.Powell, R. A. Movements, home ranges, activity, and dispersal. In Carnivore Ecology and Conservation (eds Boitani, L. & Powell, R. A.) 188–217 (Oxford University Press, 2012).Chapter 

Google Scholar 
41.Figueiredo, A. et al. Reintroduction of the Iberian lynx (Lynx pardinus): A preliminary case study in Extremadura, Spain. J. Ethol. 37, 343–351 (2019).Article 

Google Scholar 
42.Sarmento, P., Carrapato, C., Eira, C. & Silva, J. P. Spatial organization and social relations in a reintroduced population of endangered Iberian Lynx Lynx pardinus. Oryx 53, 344–355 (2017).Article 

Google Scholar 
43.Sarmento, P. & Carrapato, C. The use of spatially explicit capture-recapture models for estimating Iberian lynx abundance in a newly reintroduced population. Mamm. Biol. 98, 11–16 (2019).Article 

Google Scholar 
44.Blázquez-Cabrera, S., Ciudad, C., Gastón, A., Simón, M. A. & Saura, S. Identification of strategic corridors for restoring landscape connectivity: Application to the Iberian Lynx. Anim. Conserv. 22, 210–219 (2018).Article 

Google Scholar 
45.IBERLINCE Team Life+IBERLINCE Recovery of the historical distribution for Iberian Lynx (Lynx pardinus) in Spain and Portugal. (LIFE10NAT/ES/570) [WWW Document]. http://www.iberlince.eu/index.php/eng/project/description (2018).46.GAAS—Grupo Asesor de Aspectos Sanitarios del Lince Ibérico. Manual sanitario del lince ibérico. Versión 2.1. http://www.lynxexsitu.es/ficheros/documentos_pdf/85/Manual_Sanitario_Lince_Ib_2014.pdf (2014)47.White, G. C. & Garrott, R. A. Analysis of Wildlife Radio Tracking Data (Harcourt Brace Jovanovich, 1990).
Google Scholar 
48.Laver, P. N. & Kelly, M. J. A critical review of home range studies. J. Wildl. Manag. 72, 290–298 (2008).Article 

Google Scholar 
49.Kovach, W. L. Oriana-Circular Statistics for Windows, v.4 (Kovach Computing Services, 2011).
Google Scholar 
50.Cole, L. C. The measurement of interspecific association. Ecology 30, 411–424 (1949).Article 

Google Scholar 
51.Ferreira, C., Paupério, J. & Alves, P. C. The usefulness of field data and hunting statistics in the assessment of wild rabbit (Oryctolagus cuniculus) conservation status in Portugal. Wildl. Res. 37, 223–229 (2010).Article 

Google Scholar 
52.Monterroso, P. et al. Disease-mediated bottom-up regulation: An emergent virus affects a keystone prey, and alters the dynamics of trophic webs. Sci. Rep. 6, 36072 (2016).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
53.Signer, J. & Balkenhol, N. Reproducible home ranges (rhr): A new, user-friendly R package for analyses of wildlife telemetry data. Wildl. Soc. Bull. 39, 358–363 (2015).Article 

Google Scholar 
54.Calenge, C. The package “adehabitat” for the R software: A tool for the analysis of space and habitat use by animals. Ecol. Model. 197, 516–519 (2006).Article 

Google Scholar 
55.R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing. http://www.R-project.org/ (2020).56.Bates, D., Mächler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Soft. 67, 48 (2015).Article 

Google Scholar 
57.Barton, K. MuMIn—Multimodel Inference. R package version 1.43.6. http://CRAN.R-Project.org/package=MuMin (2019).58.Weise, F. J. et al. Cheetahs (Acinonyx jubatus) running the gauntlet: An evaluation of translocations into free-range environments in Namibia. PeerJ 3, e1346 (2015).PubMed 
PubMed Central 
Article 

Google Scholar 
59.Briers-Louw, W. D., Verschueren, S. & Leslie, A. J. Big cats return to Majete Wildlife Reserve, Malawi: Evaluating reintroduction success. Afr. J. Wildl. Res. 49, 34–50 (2019).
Google Scholar 
60.Yiu, S. et al. Early post-release movement of reintroduced lions (Panthera leo) in Dinokeng Game Reserve, Gauteng, South Africa. Eur. J. Wildl. Res. 61, 861–870 (2015).Article 

Google Scholar 
61.Armstrong, D. & Seddon, P. Directions in reintroduction biology. Trends Ecol. Evol. 23, 20–25 (2008).PubMed 
Article 
PubMed Central 

Google Scholar 
62.Jule, K. R., Leaver, L. A. & Lea, S. E. G. The effects of captive experience on reintroduction survival in carnivores: A review and analysis. Biol. Conserv. 141, 355–363 (2008).Article 

Google Scholar 
63.Reeves, J., Smith, C., Dierenfeld, E. S. & Whitehouse-Tedd, K. Captivity-induced metabolic programming in an endangered felid: Implications for species conservation. Sci. Rep. 10, 3630 (2020).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
64.Jackson, C. R., Groom, R. J., Jordan, N. R. & McNutt, J. W. The effect of relatedness and pack size on territory overlap in African wild dogs. Mov. Ecol. 5, 10 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
65.Jiménez, J. et al. Restoring apex predators can reduce mesopredator abundances. Biol. Conserv. 238, 108234 (2019).Article 

Google Scholar 
66.Schaub, M., Zink, R., Beissmann, H., Sarrazin, F. & Arlettaz, R. When to end releases in reintroduction programmes: Demographic rates and population viability analysis of Bearded Vultures in the Alps. J. Appl. Ecol. 46, 92–100 (2009).Article 

Google Scholar 
67.Bertolero, A., Pretus, J. L. & Oro, D. The importance of including survival release when assessing viability in reptile translocations. Biol. Conserv. 217, 311–320 (2018).Article 

Google Scholar 
68.Potts, J. R., Harris, S. & Giuggioli, L. Territorial dynamics and stable home range formation for central place foragers. PLoS One 7, e34033 (2012).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
69.Gil-Sánchez, J. M. et al. The use of camera trapping for estimating Iberian Lynx (Lynx pardinus) home ranges. Eur. J. Wildl. Res. 57, 1203–1211 (2011).Article 

Google Scholar 
70.Spiegel, O., Leu, S. T., Bull, C. M. & Sih, A. What’s your move? Movement as a link between personality and spatial dynamics in animal populations. Ecol. Lett. 20, 3–18 (2017).ADS 
PubMed 
Article 
PubMed Central 

Google Scholar 
71.Bremner-Harrison, S., Prodohl, P. A. & Elwood, R. W. Behavioural trait assessment as a release criterion: Boldness predicts early death in a reintroduction programme of captive bred swift fox (Vulpes velox). Anim. Conserv. 7, 313–320 (2004).Article 

Google Scholar 
72.Devineau, O., Shenk, T. M., Doherty, P. F., White, G. C. & Kahn, R. H. Assessing release protocols for Canada Lynx reintroduction in Colorado. J. Wildl. Manag. 75, 623–630 (2011).Article 

Google Scholar 
73.Reading, R. P., Miller, B. & Shepherdson, D. The value of enrichment to reintroduction success. Zoo Biol. 32, 332–341 (2013).PubMed 
Article 
PubMed Central 

Google Scholar 
74.Hardman, B. & Moro, D. Optimising reintroduction success by delayed dispersal: Is the release protocol important for hare-wallabies?. Biol. Conserv. 128, 403–411 (2006).Article 

Google Scholar 
75.Sarkar, M. et al. Movement and home range characteristics of reintroduced tiger (Panthera tigris) population in Panna Tiger Reserve, central India. Eur. J. Wildl. Res. 62, 537–547 (2016).Article 

Google Scholar  More

Body and brain size databaseThe fossil dataset consists of the hitherto largest collection of body (n = 204) and brain size estimates (n = 166) from Homo in the past ~1.0 Ma (Fig. 1). The data on hominin body size estimates are derived from our own previous study6 plus additional estimates60 and updated chronometric ages from more recent literature. Individual body size estimates are provided by specimen in Supplementary Data 1 with data sources. The bulk of data on hominin brain sizes (endocranial volume, in cm3) is derived from recent meta-analyses7,12,13,61,62 and updated chronometric information. Specific sources of these data are indicated in Supplementary Data 2, with some assessments bearing larger errors due to the incomplete state of the crania on which they are based (e.g., Arago 21, Vértesszőlős, Zuttiyeh). Each body and brain size estimate is associated with information on estimated chronometric age (dating method and data source), geographical location (longitude and latitude), and taxonomic attribution. For the exact locations per specimen, see interactive map in Supplementary Note 1. We divided the dataset into three taxonomic units: Pleistocene H. sapiens, Neanderthals, and Mid-Pleistocene Homo. Whereas hypodigms of H. sapiens and Neanderthal remains are generally agreed upon, we use “Mid-Pleistocene Homo” as a strictly analytical unit to denote African and European Middle Pleistocene hominins that predate Neanderthals and are not assigned to Homo naledi, between ~800 and 130 ka. We refrained from further division of this group due to the often fragmentary nature of fossils, unclear alpha taxonomy, and small sample size. Analyses performed within these taxonomic units minimize phylogenetic effects of, e.g., significantly different brain sizes (e.g., ref. 2). Specimens from H. naledi and Homo floresiensis had to be excluded from this analysis as for each taxon they derive from a single location and age bracket, precluding assessment of paleoclimatic variation. Limitations to the fossil datasets (see e.g., refs. 2,3,4,6) include imprecision of brain and body size estimates due to methodical and taphonomic problems, uncertainties of absolute ages that translate into uncertainties of the associated climate, and unequal sampling of hominin fossils across time and space. These limitations were incorporated into the construction of the synthetic dataset to assess the extent of their effects on the overall results for the actual fossil dataset. For all further analyses, brain and body size values were log-transformed as they increase multiplicatively.Climate reconstructionsEach body and brain size estimate required corresponding estimates of relevant climatic variables. Our climate records are numerical model estimates based on global climate reconstructions for the past 1 Ma using the global climate model emulator GCMET27. The main idea behind GCMET is that global climate model (GCM) simulations of the past 120,000 years contain sufficient information about long-term climatic changes on time scales of ≥1000 years. Given that we know the external boundary conditions, we can reconstruct previous glacial–interglacial climatic changes. The Quaternary climate is largely determined by dynamics of the Northern Hemisphere ice sheets, which, in turn, are affected by orbital variations of the Earth around the Sun and variations of atmospheric CO2. Using these factors as external boundary conditions, GCMET can emulate the climate of the Quaternary in a similar way as a state-of-the-art GCM27.The atmospheric CO2 record of the past 1 Ma that we use in this study is a composite of the EPICA CO2 record from an Antarctic ice core63 and of output from a carbon cycle model (CYCLOPS)64. The EPICA record covers the past ~800 ka, whereas we use the CYCLOPS model output to cover the time up until 1.0 Ma. Orbital variations are based on calculations by Berger and Loutre65. Ice-sheet extents for the past 800 ka are based on numerical ice-sheet model output66. For the period before 800 ka, we assumed present-day ice-sheet configurations. This is an appropriate assumption given that all but one specimen of the fossil record before 800 ka in our datasets are within Africa or southeast Asia and thus far away from ice-sheet margins, with the local GCMET climate reconstructions not affected by this simplification.For each fossil site location from the body and brain size database, we extracted a time series of the relevant climate variables, see Table 1 (also Supplementary Figs. 10 and 11). The time series were used to attach the value of each climate variable to the fossil record, both for the actual fossil data as well as for the synthetic fossil datasets.Linear modelsThe null and two alternative linear models used throughout this manuscript are defined as follows. The null model simply estimates the mean for each taxonomic group, and we refer to this model as LM-T (linear model with taxa):$$Y={beta }_{0}+{beta }_{1}times {rm{taxon}}$$
(1)
Here Y corresponds to either body or brain size (or the log-transformed thereof), whereas β0 is the intercept, which is equivalent to the mean size of the reference taxon, and β1 is a factor that reflects the deviation from this mean size for a taxonomic group (thus giving independent intercepts for Mid-Pleistocene Homo, Neanderthals, or Pleistocene H. sapiens).The first alternative model contains the effect of the climate variable X (across all taxa):$$Y={beta }_{0}+{beta }_{1}times {rm{taxon}}+{beta }_{2}times X$$
(2)
Here β0 and β1 are the intercept terms, giving taxon-specific values, and β2 is the slope, which is the same across all taxa. We refer to this model as LM-TC (linear model with taxonomic differences plus a climate effect).The second alternative model takes taxonomic differences for the slope of the climate effect into account. This is done via an interaction term, β3, which acts as a modifier for the slopes (i.e., different intercepts, given by β0 and β1, and slopes, given by β2 and β3, for each taxonomic group):$$Y={beta }_{0}+{beta }_{1}times {rm{taxon}}+{beta }_{2}times X+{beta }_{3}times {rm{taxon}}times X$$
(3)
We refer to this model as LM-T*C (linear model with taxonomic differences plus a taxon-specific climate effect). The slopes, β2 and β3 in Eqs. (2) and (3), respectively, are presented in the main text in Tables 2 and 3 (for the log-transformed sizes) and in Supplementary Tables 1 and 2 (for the natural units of the sizes).Synthetic datasets and power analysisApart from determining the smallest sample size suitable to detect the effect of a given test at the desired level of significance, power analysis can also be used as a formal way to test whether a relationship between dependent and independent variables can be detected with the available data and proposed methods (i.e., linear models in our case) assuming that such a relationship exists. Before testing for any true association between local climate and the fossil record, we use such a power analysis to assess our power to detect relationships of different effect sizes given the uncertainties, for example, in body/brain sizes, dating, and climate reconstructions. We generated 1000 synthetic datasets for each of the ten climate variable associations (MAP, MAT, NPP, mean temperature of coldest quarter, mean precipitation of driest quarter, and the logarithm of their running standard deviation over a 10,000-year window) with body and brain size. For each association, we assumed a strong, a medium, and a weak relationship between size and climate.By strong, we refer to 1/4 of the maximum possible slope given by the range of the climate and size. Subsequently, medium is half the slope of the strong relationship, (1/8 maximum possible slope), and weak is half the slope of the medium relationship (1/16 maximum possible slope). For example, the strong association between MAT and body size (Bergmann’s rule) is −0.34 kg/°C, based on the above defined rule. This is close to the estimated association between temperature and body size of about −0.4 kg/°C found for modern humans in a recent study (ref. 31, their Fig. 5A). Unfortunately, there are no empirical data about other climatic relationships and body (or brain) size. For simplicity, we therefore applied the same rule of strong, medium, and weak associations for all other climate variables and for brain size.Before generating a synthetic dataset, we estimated the intercepts β0 and β1 and the slope β2 for the LM-TC model, Eq. (2). However, for the real fossil analysis we used the model with the interaction term, LM-T*C, Eq. (3). First, we looked up the climate record for each location and time from the climate time series and attached it to the respective empirical fossil records. We calculate the maximum slope from the X and Y ranges as β1 = range(Y)/range(X). Assigning an actual relationship factor, e.g., strong (=1/4), the intercept β0 can be calculated using the X- and Y-midpoints, β0 = Ymidpoint − 1/4β1Xmidpoint.For the synthetic fossil datasets, we assume an age uncertainty range of 10% (±5%) for radiocarbon-dated fossils, i.e., younger than 50 ka cal BP (e.g. ref. 67), and 20% (±10%) for fossils older than 50 ka coming from other dating methods with higher uncertainty such as luminescence, U-series, or ESR (e.g., ref. 68). Furthermore, we assume a standard error of 2 K for mean annual and mean temperature of coldest quarter. For all other climate variables, we assume a 20% error range (±10%). The 2 K and the 20% are in line with climate model biases as estimated in a recent study69. Within a taxonomic unit of the genus Homo, we assume a coefficient of variation (CV) of 7% for body size (average of intrapopulation means of 19 global Holocene hunter-gatherer populations, n = 510, data from JTS) and 3.5% for brain size (from ref. 28 populations, dataset: http://volweb.utk.edu/~auerbach/HOWL.htm; ref. 70, see also ref. 32). Previous research has demonstrated that the range of body size variation in Holocene human populations is larger than any taxonomic unit of earlier hominin and encompasses the range of variation found within earlier hominins6 and that sexual dimorphism in size among Mid-Pleistocene hominins is comparable to that of modern humans71. While there are significant differences in brain shape through recent hominin evolution, the range of size variation within Pleistocene hominin taxa remains comparable to that observed among modern humans60. These observations suggest that modeling the intrapopulation variation among hominin taxa upon modern human coefficients of variation provides a reasonable estimate of variation within hominin taxa that are often presented only by much smaller sample sizes. To create a synthetic dataset that has a mean and a variance as close to the fossil dataset, we introduced taxonomic size differences (β1 in Eq. (2)) that is based on the taxonomic differences in the mean size. This difference was estimated directly from the fossil dataset.The procedure to generate a single synthetic dataset is as follows. First, we selected a relationship strength, e.g., strong, and calculated the slope and intercepts. For each synthetic data point, we:

1.

Looked up the age and added a randomly sampled error (±5% or ±10%).

2.

Looked up the fossil site and selected the climate record from the previously calculated time series for that location and sampled age.

3.

Added a randomly sampled error (i.e., S.D. of 2 K or 20%) to the climate record. This is now the X value.

4.

Multiplied X with the slope β2 and added the intercept β1 with the respective taxonomic correction. This translates the climate record X into a size estimate Y.

5.

Added a random term to Y based on the CV, i.e., 3.5% for brain and 7% for body size.

6.

Repeated steps (1)–(5) for each fossil record and saved all locations, ages, Xs, and Ys to a file. This is a single synthetic dataset in the same format as the original fossil dataset.

We repeated this N times to generate N synthetic datasets and repeated the same procedure for the other relationship strengths, i.e., medium, and weak and for all other climate variables. Panels of exemplary synthetic datasets for body and brain sizes in comparison with the original data are shown in Supplementary Figs. 10 and 11.We use the same thinning approach as described in the main text (n = 1000) for the synthetic datasets. These are then used for a power analysis to test whether the linear relationship between any climate variable and body or brain size can be detected. We fitted both the LM-TC model, Eq. (2) (in which the slope defining the relationship between climate and size is the same for the three taxonomic groups, which can differ in their intercept), and the LM-T*C model, Eq. (3) (different slopes and intercepts for the three groups). A climate effect was deemed present if the null model had a higher AIC value compared to either of the alternative models, LM-TC or LM-T*C, (ΔAIC  > 2), i.e., LM-T, Eq. (1), in which the three groups differ in size but there is no effect of climate. Figure 2 in the main text shows the power to detect a true relationship between size and climate. Individual records are color-coded according to the AIC difference between the LM-T and the alternative models, LM-TC or LM-T*C, ranging from −2 (red) to +15 (blue) with 2 as midpoint (white).All statistical tests were undertaken in Python version 3.8.5 using the following Python packages: statsmodels 0.12 (for linear models), pandas 1.1.3 (for dataframes, reading/writing CSV/Excel files), netCDF4 1.5.3 (reading NetCDF files), matplotlib 3.3.2 (for plotting), and numpy 1.19.2 (numerics).Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

Results of our research showed that lipid accumulation in sea urchin gonads follows a periodic fluctuation, in agreement with previous observations22. The analysis of Fig. 1 suggests the key role of photoperiod in triggering and then modulating fat utilization and storage mechanisms in P. lividus gonads, while the effect of temperature in gametogenesis and spawning in echinoderms still remains uncertain5,22,28. In fact, a change in photoperiod anticipated the corresponding change in gonad total lipids content in both habitats, while the role of temperature was not very clear, since lipid changes seemed not to be associated with changes in temperature. Most likely, the combined effect of both parameters regulates reproductive cycle of sea urchins. Similar periodical trends in total lipid content were also observed in recent studies on P. lividus gonads19,20 collected in different geographical areas. For example, Rocha et al.20 reported that gonadal lipid content are likely influenced by the environmental conditions characterizing the harvest site in the Praia Norte (Portugal). In this work and in the abovementioned studies, total lipid content in gonads changed as a function of gametogenic cycle, i.e. increased until the recovery/growing stage (I–II) and then progressively decreased until the premature/mature stage (III–IV)27. In another detailed characterization of Arbacia dufresnii, Dìaz de Vivar et al.29 observed a marked dependence of the total lipid content with gonad maturation, with a significant decrease in lipid content in spawned compared to intact gonads, especially in female sea urchins.Assuming P. oceanica and H. scoparia were the main dietary sources of lipids in our study, gonad lipid content was relatively independent from dietary lipid intake, in agreement with data from other authors19,20,22. Indeed, total lipids in P. oceanica and H. scoparia were very low (approximately 1% D.W.) and seasonal variations of lipid levels in these main dietary substrates were definitely negligible. These results are further supported by the literature30,31. For a better understanding of the comparison between different scientific reports19,20,22, it should be recalled here that the displacement of periodical gametogenic cycles is strongly influenced by several environmental factors and is, therefore, dependent on the growing habitat16,32.As far as the commercial value of sea urchin gonads is considered, several reports20,33 suggest that the best harvesting period is when gonads are in the growing stage, when nutrient contents (i.e. proteins, lipids and carbohydrates) are at their highest levels, and when sensorial characteristics are optimal. In fact, gonad maturation decreases the overall quality of roe, and make them more bitter and less pleasant33,34. However, it is striking that very often the official regulation on the harvest of P. lividus in Sardinia allowed collection of sea urchins in the period from November to April, when products are nutrient-poor and in the late stages of gametogenesis (i.e. pre-mature, mature and spawning stages)35.Our observations suggest that total lipids from dietary sources concentrate in the gut. The amount of lipids in these latter samples is actually always much higher than in the sea grass and macroalgae analyzed. The concentration of lipids in the gut has been already observed in other echinoderms as well36,37. These evidences suggest that digestion phenomena occurring in the gut may include the concentration of nutrients. Moreover, our data show that lipid fatty acid composition in gut is considerably consistent, regardless dietary lipid. While further studies are needed, most recent findings strongly suggest that gut flora have a role in assisting digestion and absorption of nutrients in sea urchins38. De novo synthesis of fatty acids by microbiotes, an interesting hypothesis that would especially concern the modulation of short chain fatty acids levels, should be further and specifically investigated. Based on most recent findings, a possible role of bacteria in nutrient production and processing has been postulated39. However, it should be also reckoned that other lipids may come from other dietary sources beyond the main sea grass and macroalgae (“Supplementary Material”). This latter hypothesis, however, would not explain the substantial increase observed in gut lipids, since other possible sources do not have very high lipid contents and were taken in small percentage. For example, it was previously observed in adult Strongylocentrotus intermedius that algal pellets exceeded 80–90% (wet weight) of gut contents, complemented by detritus, small animals (e.g. small crustaceans and mollusks) and non-foods (e.g. sand, shell fragments)40. Moreover, in P. lividus sampled from natural conditions in Corsica (France), 95% of the total gut content was represented by plant material41. Similarly, animal taxa in our study represented a very low percentage of the gut content, and species populating the rocky bottom, other than H. scoparia, have low lipid content and likely had little relevance on sea urchins diet. Also Murillo-Navarro and Jimenez-Guirado25, in a yearlong investigation, found that H. scoparia was the most abundant brown alga in gut contents of P. lividus.Brown algae and leaves of P. oceanica are in fact generally considered among the primary components of adult P. lividus diets1,24,25. It has been also observed by other authors that sea urchins consume all parts of P. oceanica and preferentially green leaves colonised by epiphytes1,26,42,43,44. Epiphytes were not removed from our samples before analysis.The role of gut and stomach as nutrient storage organs is generally acknowledged41,45. This is demonstrated by the almost double lipid contents found in gut than in food sources in the present investigation and by other studies36,37. As a later digestion step, lipids are selectively stored in gonads, where almost three or even four times the lipids contents found in the gut were detected. This supports the hypothesis of lipid relocation from gut to gonads, thus confirming the role of gonads as an important storage tissue for P. lividus, as was previously established by other authors22,46 and correspondingly a role in lipid metabolism can be ascribed to the digestive tract. It also further proves that the amount of fat daily introduced with diet has only a limited influence on the seasonal evolution of total lipids in gonads. Of course, nutrients and especially lipids stored in gonads serve during gametogenesis, as an energy source for developing embryos and are mobilized during pre-feeding development of larvae5. In echinoderms, indeed, nutrients provided in the eggs are needed by developing embryos and larvae.In two recent investigations on P. lividus collected along the Atlantic coast of Portugal, Rocha et al.19,20 evidenced slightly different seasonal trends. They observed both a maximum lipid content and an increase in PUFA content in gonads during the fall season. In contrast, we observed a peak in total lipids during summer, and an increase in PUFA during winter. Likely, the different climatic and environmental conditions of the Atlantic coast with respect to the Mediterranean basin (especially seawater temperatures) induce different gametogenesis cycles16, which in turn modulate the lipid balance in gonads. Gametogenic stages are in fact differently distributed along the year in ours and the cited works by Rocha et al.19,20. In general, lipid content in gonads seem to increase during the recovery (stage I) and growing (stage II) gametogenic stages27, when gonads are packed with nutritive phagocytes and only few germ cells are present.Other studies suggested that specific fatty acids found in the gonads of sea urchins may be synthesized by other tissues such as the intestine and then mobilized to the gonads47.Regardless the different food availability in the two analyzed sites, our results show a remarkable robustness of the fatty acids profile of gut contents. This is particularly interesting since they show a regulation of physiologically essential C 20:5 n-3 and C 20:4 n-6 at gut level, which seem to quite finely level out according to season, regardless the dietary contents of these fatty acids.The increase in gonad PUFA observed in both habitats during winter did not seem to correlate with substantial changes in the main taxa isolated in the gut content of the sea urchin sampled in the P. oceanica meadow, nor to relevant changes in the specimens populating the rocky bottom habitat (“Supplementary Material”). This is consistent with our previous studies21,22, which linked the phenomenon to both the cold acclimatization effect and gametogenesis. Raise in PUFA in lower temperatures allows maintaining cell membrane fluidity and, consequently, supports its functionality.The questions arise whether the lipid species contained in the food sources can be directly and selectively absorbed by sea urchin gonads and how much food habits affect gonads composition. In order to answer these questions, discussion should be directed to each relevant fatty acid.The fatty acids of glycerolipids of higher-plants chloroplasts are highly unsaturated, and the most represented fatty acid is C 18:3 (n-3)48. Instead, brown algae, such as Phaeophyceae, contain a large amount of C 20:4 (n-6) and C 20:5 (n-3)49. During our studies, the most significant difference between the fatty acid profiles of P. oceanica and H. scoparia was related to C 18:3 (n-3). According to our data, P. oceanica contained, on average, more than ten times the amount of this FA in H. scoparia.The fatty acid profile of P. oceanica described in the present study is in agreement with previous reports50,51 and confirms that lipids of P. oceanica are mainly represented by the C 18:3 (n-3), C 18:2 (n-6), and C 16:051. On the contrary, the fatty acid composition of H. scoparia seems to be quite variable considering previously published reports, although literature generally agrees on the most abundant fatty acids (i.e. C 16:0, C 18:2 n-6, C 20:5 n-3 and C 20:4 n-6)31,52.Both in rocky bottom and in P. oceanica meadows, gonadal C 18:3 (n-3) decreased when sea urchins metabolism is mainly influenced by production of gametes (from November), i.e. when gonads reached premature/mature stages, as previously observed15,22. Our data showed a decrease of C 18:3 (n-3) in gut roughly corresponding to an increase of the same FA in gonads (Fig. 4), suggesting that dietary C 18:3 (n-3) was not selectively and directly retained in gonads from the diet, but likely took active part to metabolic processes of bioconversion or is catabolized during β-oxidation of lipids.Also C 18:2 (n-6) showed a similar behaviour in our study and in other previous investigations15,20.Remarkably, C 20:5 (n-3) and C 20:4 (n-6) were the most abundant LC-PUFA in both gut and gonads, in contrast with the composition of the main dietary sources of lipids in the two habitats. In fact, while high percentages of these fatty acids were found in the brown algae H. scoparia, they were present only in very low percentages in the P. oceanica samples. In sea urchins, the fatty acid profile of diet is often scarcely reflected in gut contents and gonads53. From July to March we detected higher percentages of C 20:5 (n-3) in gonad samples collected from P. oceanica meadow than in the corresponding samples from rocky bottom. Moreover, our data clearly show that the C 20:5 (n-3) contained in either gonads and gut does not reflect seasonal variations of this FA in the main sea grass and macroalgae populating the two sites. This result supports earlier observations5,21.Beyond P. oceanica, green algae, especially C. cylindracea, represented additional dietary sources of C 20:5 (n-3) in the P. oceanica meadow. P. lividus usually feeds on brown algae and only less frequently on green algae1,15. In fact, green algae represented less than 5% of the gut content in P. oceanica meadow all year long but from October to December, when they increased from 10 to 25%. In this period, C 20:4 (n-6) and C 20:5 (n-3) in gonads reached their lowest values, but the C 20:5 (n-3) content in gut noticeably increased. After January, when sea urchin reduced feeding in green algae and again less than 5% of green algae was found in the gut content, C 20:5 (n-3) and C 20:4 (n-6) content in sea urchins gut started increasing. To explain this observation, we recall that it was found in S. droebachiensis that dietary FA were not incorporated in sea urchin tissues after short feeding experiments54, but longer experiments allowed to observe diet-related modifications in tissues36. Therefore, it is reasonable to think that nutrients are transferred from gut to gonads. Among other dietary sources of lipids, brown algae in P. oceanica meadow likely did not significantly contribute to increase LC PUFA in gut contents and gonads prior to gametogenesis, being brown algae intake almost always low in the present study.The observed increase of C 20:4 (n-6) in gonads in December was less correlated to the dietary availability of this FA, but was likely associated to cold adaptation and to the growth and maturation of gametes21. In fact, even when the main dietary source of lipids, P. oceanica, was almost completely devoid of this FA, the percentage of C 20:4 (n-6) in gonads was 10–15% and not significant increase of this FA was observed in gut contents from October to December.As for most aquatic consumers, C 20:5 (n-3) and C 20:4 (n-6) can be selectively retained in gonads from dietary sources or accumulated through the conversion of other essential 18-carbon FA.Since we found similar amount of C 20:5 (n-3) C 20:5 (n-3) and C 20:4 (n-6) in P. lividus gonads and gut contents and these values were much higher than in dietary sources, retention or biosynthesis should have occurred already at intestinal level, as previously suggested for other echinoderms36,37,47. As previously hypothesized for Strongylocentrotus intermedius, likely these FA were transferred to gonads after being processed and stored in the digestive tract47.Recently, Kabeya et al.55 found that P. lividus possesses desaturases that are able to convert C 18:3 (n-3) and C 18:2 (n-6) into C 20:5 (n-3) and C 20:4 (n-6), respectively. Han et al.47 characterized the expression of fatty acid desaturases (SiFad1) in different tissues of S. intermedius and concluded that the highest expression is in the intestine, while gonads have lower expression level. Therefore, while retention from diet and biosynthesis from C18 precursors of essential lipid species such as C 20:5 (n-3) and C 20:4 (n-6) might occur already in the gut36,37,41,45, also gonads might possess some, likely lower, biosynthetic functions. Kabeya et al.55 did not specifically quantify the expression of desaturases in different tissues of P. lividus, therefore further research in this sense would be beneficial.It should be mentioned that sex-induced difference of fatty acid profiles of sea urchin gonads were not studied in the present work, but males and females specimens were pooled together. Some previous reports have evidenced differences in lipid classes and fatty acids profiles between sexes15,29, while other studies did not spot statistically significant gender-related discrepancies5. Fatty acids profiles of gonads are likely to be related by sea urchin gender, but it is reasonable to believe that such differences would not disprove the aforementioned considerations on lipid storage and metabolism at gut and at gonad level. In particular, the differences in C 18:3 n-3, C 18:2 n-6, C 20:4 n-6 and C 20:5 n-3 found in previous studies between male and female gonads were quite low (maximum 2–4% of total FAME). Gender differences are ascribable to the increasing presence of lipid-rich gametes (oocytes or sperm) during the gonad maturation period. Also differences in lipid classes are expected in this period, being triglycerides mainly present in female gametes29,56. According to previous reports, during the reproductive period females of both P. lividus and Arbacia lixula showed lower proportions of 20:4n-6, while 20:5n-3 was higher in males of P. lividus and in females of A. lixula56. In P. lividus, such differences were found to be very limited for 20:4n-6 and 20:5n-3 (0.1% and 1.3%, respectively, between mean values of total FAME percentage)56. Also in Arbacia dufresnii the differences between male and female intact gonads for 20:4n-6, while 20:5n-3 were found to be not very important, but both fatty acids seemed to be slightly more concentrated in male tissues29.In any case, the present study confirms that during maturation stages of gonads, when their nutritive content decreases20,22, C 20:5 (n-3) and C 20:4 (n-6) levels increase, and so does their nutritional quality. C 20:5 (n-3) consumption is in fact associated to reduced risk of several chronic diseases57. At the same time, previous reports showed that the best commercial value of sea urchin gonads is before the onset of gametogenesis20,33. These results are quite relevant not only because they allow to deepen the knowledge of the metabolic response of sea urchin P. lividus to season and diet, but also for both improving echinoculture practices and guiding relevant policies directed to regulate the harvest of wild populations. Changes in the concentration of biochemical components in the gonads of sea urchins impact their sensory quality20,33,34. In particular, gonads in their mature stages were described as more bitter34 and of lower quality overall33 than when they are in the growing stage. On the other hand, gonads in the growing stage reach the highest contents of nutrients (protein, fat, carbohydrates)20. Harvest of wild sea urchin during the reproductive time should be avoided, and this is particularly important for an endangered species such as P. lividus. Echinoculture could provide sea urchin roe for which the harvest time should be carefully scheduled as a function of analytical quality parameters and based on expected use.In conclusion, P. oceanica and H. scoparia, primarily constituted P. lividus diet in two contrasting sites within the same geographical area. Green algae, especially C. cylindracea, supplemented sea urchin diet in the P. oceanica meadow prior to gametogenesis, demonstrating the ability of P. lividus to select their diet according to requirements. Total lipid content in gonads changed periodically as a function of gametogenic cycle, being relatively independent from dietary lipid intake and showing a maximum during the growing stage and a minimum in mature gonads. Fatty acid profiles of P. oceanica and H. scoparia were significantly different from each other throughout the year. C 18:3 (n-3) was the main differential dietary marker in P. lividus gonads and gut contents. The main PUFA of P. lividus gonads, C 20:5 (n-3) and C 20:4 (n-6) were associated to increased consumption of green algae in P. oceanica meadow. LC-PUFA were selectively allocated in gonads as a function of reproductive cycle. Conversion of C 18:3 (n-3) to C 20:5 (n-3) and of C 18:2 (n-6) to C 20:4 (n-6) at gut level cannot be excluded, although further research in this sense is desirable. It is worth to note that harvest is generally allowed in Sardinia during gonad maturation, when main nutrients (lipids, carbohydrates, proteins) are at lowest level and also the sensory quality of roe is low, but gonads are rich in healthy LC-PUFA. Our results suggest that rearing of P. lividus would be possible with diets very poor in LC-PUFA given a supplement of this nutrients is provided prior to gametogenesis, when gonads are in the growing/premature stages. More

1.Robbins, C. T. & Cunha, T. J. Wildlife Feeding and Nutrition (Elsevier Science, 2014).
Google Scholar 
2.Murray, M. H., Becker, D. J., Hall, R. J. & Hernandez, S. M. Wildlife health and supplemental feeding: A review and management recommendations. Biol. Conserv. 204, 163–174 (2016).Article 

Google Scholar 
3.Barboza, P. S., Parker, K. L., & Hume, I. D. Integrative Wildlife Nutrition (Springer, 2009).Book 

Google Scholar 
4.Nyhus, P. J. Human-wildlife conflict and coexistence. Annu. Rev. Environ. Resour. 41, 143–171 (2016).Article 

Google Scholar 
5.Baynham-Herd, Z., Redpath, S., Bunnefeld, N. & Keane, A. Predicting intervention priorities for wildlife conflicts. Conserv. Biol. 34, 232–243 (2020).PubMed 
Article 

Google Scholar 
6.Treves, A. & Santiago-Ávila, F. J. Myths and assumptions about human-wildlife conflict and coexistence. Conserv. Biol. 34, 811–818 (2020).PubMed 
Article 
PubMed Central 

Google Scholar 
7.Bojarska, K. & Selva, N. Spatial patterns in brown bear Ursus arctos diet: The role of geographical and environmental factors: Biogeographical variation in brown bear diet. Mammal Rev. 42, 120–143 (2012).Article 

Google Scholar 
8.Kavčič, I. et al. Fast food bears: Brown bear diet in a human-dominated landscape with intensive supplemental feeding. Wildl. Biol. 21, 1–8 (2015).Article 

Google Scholar 
9.Cozzi, G. et al. Anthropogenic food resources foster the coexistence of distinct life history strategies: Year-round sedentary and migratory brown bears. J. Zool. 300, 142–150 (2016).Article 

Google Scholar 
10.Lewis, D. L. et al. Foraging ecology of black bears in urban environments: Guidance for human-bear conflict mitigation. Ecosphere 6, art141 (2015).ADS 
Article 

Google Scholar 
11.Naves, J., Fernández-Gil, A., Rodríguez, C. & Delibes, M. Brown bear food habits at the border of its range: A long-term study. J. Mammal. 87, 899–908 (2006).Article 

Google Scholar 
12.Rodríguez, C., Naves, J., Fernández-Gil, A., Obeso, J. R. & Delibes, M. Long-term trends in food habits of a relict brown bear population in northern Spain: The influence of climate and local factors. Environ. Conserv. 34, 36–44 (2007).Article 

Google Scholar 
13.Ciucci, P., Tosoni, E., Di Domenico, G., Quattrociocchi, F. & Boitani, L. Seasonal and annual variation in the food habits of Apennine brown bears, central Italy. J. Mammal. 95, 572–586 (2014).Article 

Google Scholar 
14.Reynolds-Hogland, M. J., Pacifici, L. B. & Mitchell, M. S. Linking resources with demography to understand resource limitation for bears: Linking resources and demography. J. Appl. Ecol. 44, 1166–1175 (2007).Article 

Google Scholar 
15.Robbins, C. T., Schwartz, C. C. & Felicetti, L. A. Nutritional ecology of ursids: A review of newer methods and management implications. Ursus 15, 161–171 (2004).Article 

Google Scholar 
16.Can, Ö. E., D’Cruze, N., Garshelis, D. L., Beecham, J. & Macdonald, D. W. Resolving human-bear conflict: A global survey of countries, experts, and key factors: Human-bear conflict. Conserv. Lett. 7, 501–513 (2014).Article 

Google Scholar 
17.Hobson, K. A., McLellan, B. N. & Woods, J. G. Using stable carbon (δ 13C) and nitrogen (δ 15N) isotopes to infer trophic relationships among black and grizzly bears in the upper Columbia River basin, British Columbia. Can. J. Zool. 78, 1332–1339 (2000).Article 

Google Scholar 
18.Mowat, G. & Heard, D. C. Major components of grizzly bear diet across North America. Can. J. Zool. 84, 473–489 (2006).CAS 
Article 

Google Scholar 
19.Ben-David, M., Titus, K. & Beier, L. R. Consumption of salmon by Alaskan brown bears: A trade-off between nutritional requirements and the risk of infanticide?. Oecologia 138, 465–474 (2004).ADS 
PubMed 
Article 
PubMed Central 

Google Scholar 
20.Hopkins, J. B. et al. Stable isotopes to detect food-conditioned bears and to evaluate human-bear management. J. Wildl. Manag. 76, 703–713 (2012).Article 

Google Scholar 
21.Hata, A. et al. Stable isotope and DNA analyses reveal the spatial distribution of crop-foraging brown bears. J. Zool. 303, 207–217 (2017).Article 

Google Scholar 
22.Hilderbrand, G. V., Jenkins, S. G., Schwartz, C. C., Hanley, T. A. & Robbins, C. T. Effect of seasonal differences in dietary meat intake on changes in body mass and composition in wild and captive brown bears. Can. J. Zool. 77, 1623–1630 (1999).Article 

Google Scholar 
23.Rode, K. D., Farley, S. D. & Robbins, C. T. Sexual dimorphism, reproductive strategy, and human activities determine resource use by brown bears. Ecology 87, 2636–2646 (2006).PubMed 
Article 

Google Scholar 
24.Hilderbrand, G. V. et al. Use of stable isotopes to determine diets of living and extinct bears. Can. J. Zool. 74, 2080–2088 (1996).Article 

Google Scholar 
25.Murray, M. H., Fassina, S., Hopkins, J. B., Whittington, J. & St. Clair, C. C. Seasonal and individual variation in the use of rail-associated food attractants by grizzly bears (Ursus arctos) in a national park. PLoS ONE 12, e0175658 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
26.Mizukami, R. N., Goto, M., Izumiyama, S., Hayashi, H. & Yoh, M. Estimation of feeding history by measuring carbon and nitrogen stable isotope ratios in hair of Asiatic black bears. Ursus 16, 93–101 (2005).Article 

Google Scholar 
27.Mizukami, R. N. et al. Temporal diet changes recorded by stable isotopes in Asiatic black bear (Ursus thibetanus) hair. Isotopes Environ. Health Stud. 41, 87–94 (2005).CAS 
PubMed 
Article 

Google Scholar 
28.Hopkins, J. B. & Kurle, C. M. Measuring the realized niches of animals using stable isotopes: From rats to bears. Methods Ecol. Evol. 7, 210–221 (2016).Article 

Google Scholar 
29.Layman, C. A. et al. Applying stable isotopes to examine food-web structure: An overview of analytical tools. Biol. Rev. 87, 545–562 (2012).PubMed 
Article 

Google Scholar 
30.Careddu, G., Calizza, E., Costantini, M. L. & Rossi, L. Isotopic determination of the trophic ecology of a ubiquitous key species—The crab Liocarcinus depurator (Brachyura: Portunidae). Estuar. Coast. Shelf Sci. 191, 106–114 (2017).ADS 
CAS 
Article 

Google Scholar 
31.Blasi, M. F. et al. Assessing resource use patterns of Mediterranean loggerhead sea turtles Caretta caretta (Linnaeus, 1758) through stable isotope analysis. Eur. Zool. J. 85, 71–87 (2018).Article 
CAS 

Google Scholar 
32.Cicala, D. et al. Spatial variation in the feeding strategies of Mediterranean fish: Flatfish and mullet in the Gulf of Gaeta (Italy). Aquat. Ecol. 53, 529–541 (2019).CAS 
Article 

Google Scholar 
33.Bearhop, S., Adams, C. E., Waldron, S., Fuller, R. A. & Macleod, H. Determining trophic niche width: A novel approach using stable isotope analysis: Stable isotopes as measures of niche width. J. Anim. Ecol. 73, 1007–1012 (2004).Article 

Google Scholar 
34.Newsome, S. D., Martinez del Rio, C., Bearhop, S. & Phillips, D. L. A niche for isotopic ecology. Front. Ecol. Environ. 5, 429–436 (2007).Article 

Google Scholar 
35.Hopkins, J. B. & Ferguson, J. M. Estimating the diets of animals using stable isotopes and a comprehensive Bayesian mixing model. PLoS ONE 7, e28478 (2012).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
36.Phillips, D. L. Converting isotope values to diet composition: The use of mixing models. J. Mammal. 93, 342–352 (2012).Article 

Google Scholar 
37.Madeira, F. et al. Stable carbon and nitrogen isotope signatures to determine predator dispersal between alfalfa and maize. Biol. Control 77, 66–75 (2014).Article 

Google Scholar 
38.García-Vázquez, A., Pinto-Llona, A. C. & Grandal-d’Anglade, A. Brown bear (Ursus arctos L.) palaeoecology and diet in the Late Pleistocene and Holocene of the NW of the Iberian Peninsula: A study on stable isotopes. Quat. Int. 481, 42–51 (2018).Article 

Google Scholar 
39.Hilderbrand, G. V. et al. The importance of meat, particularly salmon, to body size, population productivity, and conservation of North American brown bears. Can. J. Zool. 77, 132–138 (1999).Article 

Google Scholar 
40.Felicetti, L. A. et al. Use of sulfur and nitrogen stable isotopes to determine the importance of whitebark pine nuts to Yellowstone grizzly bears. Can. J. Zool. 81, 763–770 (2003).Article 

Google Scholar 
41.Schwartz, C. C. et al. Use of isotopic sulfur to determine whitebark pine consumption by Yellowstone bears: A reassessment. Wildl. Soc. Bull. 38, 664–670 (2014).Article 

Google Scholar 
42.Hopkins, J. B., Koch, P. L., Ferguson, J. M. & Kalinowski, S. T. The changing anthropogenic diets of American black bears over the past century in Yosemite National Park. Front. Ecol. Environ. 12, 107–114 (2014).Article 

Google Scholar 
43.Bentzen, T. W., Shideler, R. T. & O’Hara, T. M. Use of stable isotope analysis to identify food-conditioned grizzly bears on Alaska’s North Slope. Ursus 25, 14 (2014).Article 

Google Scholar 
44.Teunissen van Manen, J. L., Muller, L. I., Li, Z., Saxton, A. M. & Pelton, M. R. Using stable isotopes to assess dietary changes of American black bears from 1980 to 2001. Isotopes Environ. Health Stud. 50, 382–398 (2014).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
45.Braunstein, J. L., Clark, J. D., Williamson, R. H. & Stiver, W. H. Black bear movement and food conditioning in an exurban landscape. J. Wildl. Manag. 84, 1038–1050 (2020).Article 

Google Scholar 
46.Narita, R., Mano, T., Yokoyama, R. & Takayanagi, A. Variation in maize consumption by brown bears (Ursus arctos ) in two coastal areas of Hokkaido, Japan. Mammal Study 36, 33–39 (2011).Article 

Google Scholar 
47.Matsubayashi, J., Morimoto, J., Mano, T., Aryal, A. & Nakamura, F. Using stable isotopes to understand the feeding ecology of the Hokkaido brown bear (Ursus arctos) in Japan. Ursus 25, 87–97 (2014).Article 

Google Scholar 
48.Javornik, J. et al. Effects of ethanol storage and lipids on stable isotope values in a large mammalian omnivore. J. Mammal. 100, 150–157 (2019).Article 

Google Scholar 
49.Pauli, J. N., Whiteman, J. P., Riley, M. D. & Middleton, A. D. Defining noninvasive approaches for sampling of vertebrates. Conserv. Biol. 24, 349–352 (2010).PubMed 
Article 
PubMed Central 

Google Scholar 
50.Ueda, M. & Bell, L. S. Assessing dual hair sampling for isotopic studies of grizzly bears. Rapid Commun. Mass Spectrom. 33, 1475–1480 (2019).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
51.Inger, R. & Bearhop, S. Applications of stable isotope analyses to avian ecology: Avian stable isotope analysis. Ibis 150, 447–461 (2008).Article 

Google Scholar 
52.Lerner, J. E. et al. Evaluating the use of stable isotope analysis to infer the feeding ecology of a growing US gray seal (Halichoerus grypus) population. PLoS ONE 13, e0192241 (2018).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
53.Woods, J. G. et al. Genetic tagging of free-ranging black and brown bears. Wildl. Soc. Bull. 1973–2006(27), 616–627 (1999).
Google Scholar 
54.Ciucci, P. et al. Estimating abundance of the remnant Apennine brown bear population using multiple noninvasive genetic data sources. J. Mammal. 96, 206–220 (2015).Article 

Google Scholar 
55.Kendall, K. C. et al. Using bear rub data and spatial capture-recapture models to estimate trend in a brown bear population. Sci. Rep. 9, 16804 (2019).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
56.Kendall, K. C. et al. Grizzly bear density in glacier National Park, Montana. J. Wildl. Manag. 72, 1693–1705 (2008).Article 

Google Scholar 
57.Darimont, C. T. & Reimchen, T. E. Intra-hair stable isotope analysis implies seasonal shift to salmon in gray wolf diet. Can. J. Zool. 80, 1638–1642 (2002).Article 

Google Scholar 
58.Ayliffe, L. K. et al. Turnover of carbon isotopes in tail hair and breath CO2 of horses fed an isotopically varied diet. Oecologia 139, 11–22 (2004).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
59.Schwertl, M., Auerswald, K. & Schnyder, H. Reconstruction of the isotopic history of animal diets by hair segmental analysis. Rapid Commun. Mass Spectrom. 17, 1312–1318 (2003).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
60.Jones, E. S., Heard, D. C. & Gillingham, M. P. Temporal variation in stable carbon and nitrogen isotopes of grizzly bear guardhair and underfur. Wildl. Soc. Bull. 34, 1320–1325 (2006).Article 

Google Scholar 
61.Jacoby, M. E. et al. Trophic Relations of brown and black bears in several western North American ecosystems. J. Wildl. Manag. 63, 921 (1999).Article 

Google Scholar 
62.Jimbo, M. et al. Hair growth in brown bears and its application to ecological studies on wild bears. Mammal Study 45, 1–9 (2020).Article 

Google Scholar 
63.Mosbacher, J. B., Michelsen, A., Stelvig, M., Hendrichsen, D. K. & Schmidt, N. M. Show me your rump hair and I will tell you what you ate—the dietary history of muskoxen (Ovibos moschatus) revealed by sequential stable isotope analysis of guard hairs. PLoS ONE 11, e0152874 (2016).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
64.Hopkins, J. B., Ferguson, J. M., Tyers, D. B. & Kurle, C. M. Selecting the best stable isotope mixing model to estimate grizzly bear diets in the Greater Yellowstone Ecosystem. PLoS ONE 12, e0174903 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
65.Mowat, G., Curtis, P. J. & Lafferty, D. J. R. The influence of sulfur and hair growth on stable isotope diet estimates for grizzly bears. PLoS ONE 12, e0172194 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
66.Adams, M. S. et al. Intrapopulation diversity in isotopic niche over landscapes: Spatial patterns inform conservation of bear–salmon systems. Ecosphere 8, e01843 (2017).Article 

Google Scholar 
67.Reimchen, T. E. & Klinka, D. R. Niche differentiation between coat colour morphs in the Kermode bear (Ursidae) of coastal British Columbia. Biol. J. Linn. Soc. 122, 274–285 (2017).Article 

Google Scholar 
68.Kaczensky, P. et al. Status, Management and Distribution of Large Carnivores—Bear, Lynx, Wolf & Wolverine—in Europe (Verlag nicht ermittelbar, 2013).
Google Scholar 
69.Rondinini, C., Battistoni, A., Peronace, V. & Teofili, C. Lista Rossa IUCN dei Vertebrati Italiani. Comitato Italiano IUCN e Ministero dell’Ambiente e del Mare, Roma 56, (2013).70.Ciucci, P. & Boitani, L. The Apennine brown bear: A critical review of its status and conservation problems. Ursus 19, 130–145 (2008).Article 

Google Scholar 
71.Ciucci, P. et al. Distribution of the brown bear (Ursus arctos marsicanus) in the Central Apennines, Italy, 2005–2014. Hystrix Ital. J. Mammal. 28, 86–91 (2017).
Google Scholar 
72.Maiorano, L., Chiaverini, L., Falco, M. & Ciucci, P. Combining multi-state species distribution models, mortality estimates, and landscape connectivity to model potential species distribution for endangered species in human dominated landscapes. Biol. Conserv. 237, 19–27 (2019).Article 

Google Scholar 
73.Benazzo, A. et al. Survival and divergence in a small group: The extraordinary genomic history of the endangered Apennine brown bear stragglers. Proc. Natl. Acad. Sci. 114, E9589–E9597 (2017).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
74.Gervasi, V. & Ciucci, P. Demographic projections of the Apennine brown bear population Ursus arctos marsicanus (Mammalia: Ursidae) under alternative management scenarios. Eur. Zool. J. 85, 242–252 (2018).Article 

Google Scholar 
75.Clevenger, A. P., Purroy, F. J. & Pelton, M. R. Food habits of brown bears (Ursus arctos) in the Cantabrian Mountains, Spain. J. Mammal. 73, 415–421 (1992).Article 

Google Scholar 
76.Servheen, C. Conservation of small bear populations through strategic planning. Ursus 10, 67–73 (1998).
Google Scholar 
77.Tosoni, E., Mei, M. & Ciucci, P. Ants as food for Apennine brown bears. Eur. Zool. J. 85, 342–348 (2018).Article 

Google Scholar 
78.Pritchard, G. T. & Robbins, C. T. Digestive and metabolic efficiencies of grizzly and black bears. Can. J. Zool. 68, 1645–1651 (1990).Article 

Google Scholar 
79.Cameron, M. D. et al. Body size plasticity in North American black and brown bears. Ecosphere 11, e03235 (2020).Article 

Google Scholar 
80.Stock, B. C. et al. Analyzing mixing systems using a new generation of Bayesian tracer mixing models. PeerJ 6, e5096 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
81.Banner, K. M., Irvine, K. M. & Rodhouse, T. J. The use of Bayesian priors in Ecology: The good, the bad and the not great. Methods Ecol. Evol. 11, 882–889 (2020).Article 

Google Scholar 
82.Lemoine, N. P. Moving beyond noninformative priors: Why and how to choose weakly informative priors in Bayesian analyses. Oikos 128, 912–928 (2019).Article 

Google Scholar 
83.Franco-Trecu, V. et al. Bias in diet determination: Incorporating traditional methods in Bayesian mixing models. PLoS ONE 8, e80019 (2013).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
84.Johnson, D. L., Henderson, M. T., Anderson, D. L., Booms, T. L. & Williams, C. T. Bayesian stable isotope mixing models effectively characterize the diet of an Arctic raptor. J. Anim. Ecol. 89, 2972–2985 (2020).PubMed 
Article 

Google Scholar 
85.Swan, G. J. F. et al. Evaluating Bayesian stable isotope mixing models of wild animal diet and the effects of trophic discrimination factors and informative priors. Methods Ecol. Evol. 11, 139–149 (2020).Article 

Google Scholar 
86.Ward, E. J., Semmens, B. X. & Schindler, D. E. Including source uncertainty and prior information in the analysis of stable isotope mixing models. Environ. Sci. Technol. 44, 4645–4650 (2010).ADS 
CAS 
PubMed 
Article 

Google Scholar 
87.Keis, M., Tammeleht, E., Valdmann, H. & Saarma, U. Ants in brown bear diet, and discovery of a new ant species for Estonia from brown bear scats. Hystrix Ital. J. Mammal. 30, 0 (2019).
Google Scholar 
88.Warlick, A. et al. Using Bayesian stable isotope mixing models and generalized additive models to resolve diet changes for fish-eating killer whales Orcinus orca. Mar. Ecol. Prog. Ser. 649, 189–200 (2020).Article 

Google Scholar 
89.Derbridge, J. J. et al. Experimentally derived δ13C and δ15N discrimination factors for gray wolves and the impact of prior information in Bayesian mixing models. PLoS ONE 10, e0119940 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
90.Chiaradia, A., Forero, M. G., McInnes, J. C. & Ramírez, F. Searching for the true diet of marine predators: Incorporating Bayesian priors into stable isotope mixing models. PLoS ONE 9, e92665 (2014).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
91.Ciucci, P., Mancinelli, S., Boitani, L., Gallo, O. & Grottoli, L. Anthropogenic food subsidies hinder the ecological role of wolves: Insights for conservation of apex predators in human-modified landscapes. Glob. Ecol. Conserv. 21, e00841 (2020).Article 

Google Scholar 
92.Galluzzi, A., Donfrancesco, V., Mastrantonio, G., Sulli, C. & Ciucci, P. Cost of coexisting with a relict large carnivore population: Impact of Apennine brown bears, 2005–2015. Animals 11, 1453 (2021).PubMed 
PubMed Central 
Article 

Google Scholar 
93.Dahle, B., Sørensen, O. J., Wedul, E. H., Swenson, J. E. & Sandegren, F. The diet of brown bears Ursus arctos in central Scandinavia: Effect of access to free-ranging domestic sheep Ovis aries. Wildl. Biol. 4, 147–158 (1998).Article 

Google Scholar 
94.Persson, I.-L., Wikan, S., Swenson, J. E. & Mysterud, I. The diet of the brown bear Ursus arctos in the Pasvik Valley, northeastern Norway. Wildl. Biol. 7, 27–37 (2001).Article 

Google Scholar 
95.Welch, C. A., Keay, J., Kendall, K. C. & Robbins, C. T. Constraints on frugivory by bears. Ecology 78, 1105–1119 (1997).Article 

Google Scholar 
96.Rode, K. D., Robbins, C. T. & Shipley, L. A. Constraints on herbivory by grizzly bears. Oecologia 128, 62–71 (2001).ADS 
PubMed 
Article 

Google Scholar 
97.Robbins, C. T. et al. Optimizing protein intake as a foraging strategy to maximize mass gain in an omnivore. Oikos 116, 1675–1682 (2007).Article 

Google Scholar 
98.Orlandi, L. et al. The effects of nitrogen pollutants on the isotopic signal (δ 15N) of Ulva lactuca: Microcosm experiments. Mar. Pollut. Bull. 115, 429–435 (2017).CAS 
PubMed 
Article 

Google Scholar 
99.Fiorentino, F. et al. Epilithon δ15N signatures indicate the origins of nitrogen loading and its seasonal dynamics in a volcanic Lake. Ecol. Indic. 79, 19–27 (2017).Article 
CAS 

Google Scholar 
100.Noyce, K. V., Kannowski, P. B. & Riggs, M. R. Black bears as ant-eaters: Seasonal associations between bear myrmecophagy and ant ecology in north-central Minnesota. Can. J. Zool. 75, 1671–1686 (1997).Article 

Google Scholar 
101.Auger, J., Ogborn, G. L., Pritchett, C. L. & Black, H. L. selection of ants by the American black bear (Ursus americanos). West. North Am. Nat. 64, 166–174 (2004).
Google Scholar 
102.Fujiwara, S., Koike, S., Yamazaki, K., Kozakai, C. & Kaji, K. Direct observation of bear myrmecophagy: Relationship between bears’ feeding habits and ant phenology. Mamm. Biol. 78, 34–40 (2013).Article 

Google Scholar 
103.Elgmork, K. & Kaasa, J. Food habits and foraging of the brown bear Ursus arctos in central South Norway. Ecography 15, 101–110 (1992).Article 

Google Scholar 
104.Swenson, J. E., Jansson, A., Riig, R. & Sandegren, F. Bears and ants: Myrmecophagy by brown bears in central Scandinavia. Can. J. Zool. 77, 551–561 (1999).Article 

Google Scholar 
105.Costello, C. M. et al. Diet and macronutrient optimization in wild ursids: A comparison of grizzly bears with sympatric and allopatric black bears. PLoS ONE 11, e0153702 (2016).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
106.Stenset, N. E. et al. Seasonal and annual variation in the diet of brown bears Ursus arctos in the boreal forest of southcentral Sweden. Wildl. Biol. 22, 107–116 (2016).Article 

Google Scholar 
107.Eagle, T. C. & Pelton, M. R. Seasonal nutrition of black bears in the Great Smoky Mountains National Park. Bears Their Biol. Manag. 5, 94 (1983).Article 

Google Scholar 
108.Redford, K. H. & Dorea, J. G. The nutritional value of invertebrates with emphasis on ants and termites as food for mammals. J. Zool. 203, 385–395 (2009).Article 

Google Scholar 
109.Rode, K. D. & Robbins, C. T. Why bears consume mixed diets during fruit abundance. Can. J. Zool. 78, 1640–1645 (2000).Article 

Google Scholar 
110.Erlenbach, J. A., Rode, K. D., Raubenheimer, D. & Robbins, C. T. Macronutrient optimization and energy maximization determine diets of brown bears. J. Mammal. 95, 160–168 (2014).Article 

Google Scholar 
111.Charnov, E. L. Optimal foraging, the marginal value theorem. Theor. Popul. Biol. 9, 129–136 (1976).CAS 
PubMed 
MATH 
Article 
PubMed Central 

Google Scholar 
112.Mealey, S. P. The natural food habits of grizzly bears in Yellowstone National Park, 1973–74. Bears Biol. Manag. 4, 281 (1980).
Google Scholar 
113.Cicnjak, L., Huber, D., Roth, H. U., Ruff, R. L. & Vinovrski, Z. Food habits of brown bears in Plitvice Lakes National Park, Yugoslavia. Bears Biol. Manag. 7, 221 (1987).
Google Scholar 
114.Hamer, D. & Herrero, S. Grizzly bear food and habitat in the front ranges of Banff National Park, Alberta. Bears Biol. Manag. 7, 199 (1987).
Google Scholar 
115.McLellan, B. N. & Hovey, F. W. The diet of grizzly bears in the Flathead River drainage of southeastern British Columbia. Can. J. Zool. 73, 704–712 (1995).Article 

Google Scholar 
116.Sikes, R. S. & Gannon, W. L. Guidelines of the American Society of Mammalogists for the use of wild mammals in research. J. Mammal. 92, 235–253 (2011).Article 

Google Scholar 
117.Piovesan, G., Bernabei, M., Di Filippo, A., Romagnoli, M. & Schirone, B. A long-term tree ring beech chronology from a high-elevation old-growth forest of Central Italy. Dendrochronologia 21, 13–22 (2003).Article 

Google Scholar 
118.Mancinelli, S., Boitani, L. & Ciucci, P. Determinants of home range size and space use patterns in a protected wolf (Canis lupus) population in the central Apennines, Italy. Can. J. Zool. 96, 828–838 (2018).Article 

Google Scholar 
119.Gervasi, V. et al. Estimating survival in the Apennine brown bear accounting for uncertainty in age classification. Popul. Ecol. 59, 119–130 (2017).Article 

Google Scholar 
120.Hopkins, J. B. et al. A proposed lexicon of terms and concepts for human–bear management in North America. Ursus 21, 154–168 (2010).Article 

Google Scholar 
121.Costantini, M. L., Calizza, E. & Rossi, L. Stable isotope variation during fungal colonisation of leaf detritus in aquatic environments. Fungal Ecol. 11, 154–163 (2014).Article 

Google Scholar 
122.Rossi, L., di Lascio, A., Carlino, P., Calizza, E. & Costantini, M. L. Predator and detritivore niche width helps to explain biocomplexity of experimental detritus-based food webs in four aquatic and terrestrial ecosystems. Ecol. Complex. 23, 14–24 (2015).Article 

Google Scholar 
123.Ponsard, S. & Arditi, R. Detecting omnivory with δ15N. Trends Ecol. Evol. 16, 20–21 (2001).Article 

Google Scholar 
124.Smith, J. A., Mazumder, D., Suthers, I. M. & Taylor, M. D. To fit or not to fit: Evaluating stable isotope mixing models using simulated mixing polygons. Methods Ecol. Evol. 4, 612–618 (2013).Article 

Google Scholar 
125.Vehtari, A., Gelman, A. & Gabry, J. Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Stat. Comput. 27, 1413–1432 (2017).MathSciNet 
MATH 
Article 

Google Scholar 
126.McElreath, R. Statistical rethinking: a Bayesian course with examples in R and Stan (Taylor and Francis, CRC Press, 2020).Book 

Google Scholar 
127.Stock, B., Jackson, A., Ward, E. & Venkiteswaran, J. Brianstock/Mixsiar 3.1.9. (Zenodo, 2018) https://doi.org/10.5281/ZENODO.1209993.128.Koch, P. L. & Phillips, D. L. Incorporating concentration dependence in stable isotope mixing models: A reply to Robbins, Hilderbrand and Farley (2002). Oecologia 133, 14–18 (2002).ADS 
PubMed 
Article 
PubMed Central 

Google Scholar 
129.Phillips, D. L. & Koch, P. L. Incorporating concentration dependence in stable isotope mixing models. Oecologia 130, 114–125 (2002).ADS 
PubMed 
Article 
PubMed Central 

Google Scholar 
130.Phillips, D. L. et al. Best practices for use of stable isotope mixing models in food-web studies. Can. J. Zool. 92, 823–835 (2014).Article 

Google Scholar 
131.R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, Vienna, Austria, 2020). More

1.Schimper, A. F. W., Fisher, W. R., Groom, P. & Balfour, I. B. Plant-Geography Upon a Physiological Basis. Rev. and ed. edn (Clarendon Press, 1903).2.Grime, J. & Pierce, S. The Evolutionary Strategies that Shape Ecosystems (Wiley-Blackwell, 2012).Book 

Google Scholar 
3.McGill, B., Enquist, B., Weiher, E. & Westoby, M. Rebuilding community ecology from functional traits. Trends Ecol. Evol. 21, 178–185 (2006).PubMed 
Article 
PubMed Central 

Google Scholar 
4.Chapin Iii, F. S., Bret-Harte, M., Hobbie, S. & Zhong, H. Plant functional types as predictors of transient responses of arctic vegetation to global change. J. Veg. Sci. 7, 347 (1996).Article 

Google Scholar 
5.Grime, J. P. Plant Strategies, Vegetation Processes, and Ecosystem Properties (Wiley, 2001).
Google Scholar 
6.Lavorel, S. & Garnier, E. Aardvarck to Zyzyxia-functional groups across kingdoms. New Phytol. 149, 360–363 (2001).PubMed 
Article 
PubMed Central 

Google Scholar 
7.Guo, W. et al. The role of adaptive strategies in plant naturalization. Ecol. Lett. 21, 1380–1389 (2018).PubMed 
Article 

Google Scholar 
8.Pierce, S., Luzzaro, A., Caccianiga, M., Ceriani, R. & Cerabolini, B. Disturbance is the principal α-scale filter determining niche differentiation, coexistence and biodiversity in an alpine community. J. Ecol. 95, 698–706 (2007).Article 

Google Scholar 
9.Pinho, B., Tabarelli, M., Engelbrecht, B., Sfair, J. & Melo, F. Plant functional assembly is mediated by rainfall and soil conditions in a seasonally dry tropical forest. Basic Appl. Ecol. (2019).10.Wang, J. et al. Plant community ecological strategy assembly response to yak grazing in an alpine meadow on the eastern Tibetan Plateau. Land Degrad. Dev. 29, 2920–2931 (2018).Article 

Google Scholar 
11.Barba-Escoto, L., Ponce-Mendoza, A., García-Romero, A. & Calvillo-Medina, R. P. Plant community strategies responses to recent eruptions of Popocatépetl volcano, Mexico. J. Veg. Sci. 30, 375–385 (2019).Article 

Google Scholar 
12.Diaz, S., Cabido, M. & Casanoves, F. Plant functional traits and environmental filters at a regional scale. J. Veg. Sci. 9, 113–122 (1998).Article 

Google Scholar 
13.Kelly, R. et al. Climatic and evolutionary contexts are required to infer plant life history strategies from functional traits at a global scale. Ecol. Lett. 24, 970 (2021).PubMed 
Article 

Google Scholar 
14.Odum, E. P. The strategy of ecosystem development. Science 164, 262–270 (1969).ADS 
CAS 
PubMed 
Article 

Google Scholar 
15.Reich, P. The world-wide “fast-slow” plant economics spectrum: A traits manifesto. J. Ecol. 102, 275 (2014).Article 

Google Scholar 
16.Rosado, B. H. P. & De Mattos, E. A. On the relative importance of CSR ecological strategies and integrative traits to explain species dominance at local scales. Funct. Ecol. 31, 1969 (2017).Article 

Google Scholar 
17.Raunkiær, C. The Life Forms of Plants and Statistical Plant Geography (Oxford University Press, 1934).
Google Scholar 
18.MacArthur, R. H. & Wilson, E. O. The Theory of Island Biogeography (Princeton University Press, 1967).
Google Scholar 
19.Grime, J. P. Vegetation classification by reference to strategies. Nature 250, 26–31 (1974).ADS 
Article 

Google Scholar 
20.Grime, J. P. Evidence for the existence of three primary strategies in plants and its relevance to ecological and evolutionary theory. Am. Nat. 111, 1169–1194 (1977).Article 

Google Scholar 
21.Liao, H. et al. The role of functional strategies in global plant distribution. Ecography n/a (2020).22.Pierce, S. et al. A global method for calculating plant CSR ecological strategies applied across biomes world-wide. Funct. Ecol. 31, 444–457 (2017).Article 

Google Scholar 
23.Junker, R., Lechleitner, M., Kuppler, J. & Ohler, L.-M. Interconnectedness of the Grinnellian and Eltonian niche in regional and local plant-pollinator communities. Front. Plant Sci. 10, 1371 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
24.Yu, R., Huang, J., Xu, Y., Ding, Y. & Zang, R. Plant functional niches in forests across four climatic zones: Exploring the periodic table of niches based on plant functional traits. Front. Plant Sci. 11, 841 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
25.Westoby, M. A leaf-height-seed (LHS) plant ecology strategy scheme. Plant Soil 199, 213–227 (1998).CAS 
Article 

Google Scholar 
26.Westoby, M., Falster, D., Moles, A., Vesk, P. & Wright, I. Plant ecological strategies: Some leading dimensions of variation between species. Annu. Rev. Ecol. Syst. 33, 125–159 (2002).Article 

Google Scholar 
27.Diaz, S. et al. The global spectrum of plant form and function. Nature 529, 167–171 (2016).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
28.Pierce, S. & Cerabolini, B. Plant economics and size trait spectra are both explained by one theory. (2018).29.Grime, J. P. Plant Strategies and Vegetation Processes (Wiley, 1979).
Google Scholar 
30.Grime, J. P. A comment on Loehle’s critique of the triangular model of primary plant strategies. Ecology 69, 1618–1620 (1988).Article 

Google Scholar 
31.Grime, J. et al. Integrated screening validates primary axes of specialisation in plants. Oikos 79, 259–281 (1997).Article 

Google Scholar 
32.Hodgson, J. G., Wilson, P. J., Hunt, R., Grime, J. P. & Thompson, K. Allocating C-S-R plant functional types: A soft approach to a hard problem. Oikos 85, 282–294 (1999).Article 

Google Scholar 
33.Pierce, S. & Cerabolini, B. E. L. Allocating CSR plant functional types: The use of leaf economics and size traits to classify woody and herbaceous vascular plants. Funct. Ecol. 27, 1002–1010 (2013).Article 

Google Scholar 
34.Cerabolini, B. E. L. et al. Can CSR classification be generally applied outside Britain?. Plant Ecol. 210, 253–261 (2010).Article 

Google Scholar 
35.Shipley, B. & Li, Y. An experimental test of CSR theory using a globally calibrated ordination method. PLoS ONE 12, e0175404 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
36.Rosenfield, M. F., Müller, S. C. & Overbeck, G. E. Short gradient, but distinct plant strategies: The CSR scheme applied to subtropical forests. J. Veg. Sci. 30, 984–993 (2019).Article 

Google Scholar 
37.Pyšek, P., Sádlo, J., Mandák, B. & Jarosík, V. Czech alien flora and the historical pattern of its formation: What came first to Central Europe?. Oecologia 135, 122–130 (2003).ADS 
PubMed 
Article 

Google Scholar 
38.Lambdon, P., Lloret, F. & Hulme, P. Do alien plants on Mediterranean islands tend to invade different niches from native species?. Biol. Invasions 10, 703–716 (2008).Article 

Google Scholar 
39.Dainese, M. & Bragazza, L. Plant traits across different habitats of the Italian Alps: A comparative analysis between native and alien species. Alpine Bot. 122, 11–21 (2012).Article 

Google Scholar 
40.Alexander, J. et al. Plant invasions into mountains and alpine ecosystems: Current status and future challenges. Alpine Bot. 126, 89 (2016).Article 

Google Scholar 
41.Condit, R. Tropical Forest Census Plots: Methods and Results from Barro Colorado Island, Panama and a Comparison with Other Plots (Springer, 1998).Book 

Google Scholar 
42.Pérez-Harguindeguy, N. et al. New handbook for standardised measurement of plant functional traits worldwide. Aust. J. Bot. 61, 167–234 (2013).Article 

Google Scholar 
43.Cerabolini, B. et al. Why are many anthropogenic agroecosystems particularly species-rich?. Plant Biosyst. 150, 550–557 (2014).Article 

Google Scholar 
44.Team, R. C. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2020).
Google Scholar 
45.Ferry, N. E. H. A. M. {ggtern}: Ternary diagrams using {ggplot2}. J. Stat Softw. 87, 1–17 (2018).
Google Scholar 
46.Pinheiro, J. B. D., DebRoy, S., Sarkar, D., & R Core Team. nlme: Linear and Nonlinear Mixed Effects Models (2021).47.Kassambara, A. ggpubr: ‘ggplot2’ Based Publication Ready Plots (2020).48.Diaz, S. et al. The plant traits that drive ecosystems: Evidence from three continents. J. Veg. Sci. 15, 295–304 (2004).Article 

Google Scholar 
49.Wright, I. J. et al. The worldwide leaf economics spectrum. Nature 428, 821 (2004).ADS 
CAS 
Article 

Google Scholar 
50.Parkhurst, D. F. & Loucks, O. L. Optimal leaf size in relation to environment. J. Ecol. 60, 505–537 (1972).Article 

Google Scholar 
51.Fonseca, C., Overton, J., Collins, B. & Westoby, M. Shifts in trait-combinations along rainfall and phosphorus gradients. J. Ecol. 88, 964–977 (2001).Article 

Google Scholar 
52.Hodgson, J. et al. Is leaf dry matter content a better predictor of soil fertility than specific leaf area?. Ann. Bot. 108, 1337–1345 (2011).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
53.Han, X.-W., Fang, J. Y., Reich, P., Woodward, I. & Wang, Z. Biogeography and variability of eleven mineral elements in plant leaves across gradients of climate, soil and plant functional type in China. Ecol. Lett. 14, 788–796 (2011).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
54.Ordoñez, J. et al. A global study of relationships between leaf traits, climate and soil measures of nutrient fertility. Glob. Ecol. Biogeogr. 18, 137–149 (2009).Article 

Google Scholar 
55.Bernard-Verdier, M. et al. Community assembly along a soil depth gradient: Contrasting patterns of plant trait convergence and divergence in a Mediterranean rangeland. J. Ecol. 100, 1422–1433 (2012).Article 

Google Scholar 
56.Freschet, G. et al. Global to community scale differences in the prevalence of convergent over divergent leaf trait distributions in plant assemblagesg eb_651 755..765. Glob. Ecol. Biogeogr. 20, 755–765 (2011).Article 

Google Scholar 
57.Niinemets, Ü. Global-scale climatic controls of leaf dry mass per area, density, and thickness in trees and shrubs. Ecology 82, 453–469 (2001).Article 

Google Scholar 
58.Grime, J. P. Benefits of plant diversity to ecosystems: Immediate, filter and founder effects. J. Ecol. 86, 902–910 (1998).Article 

Google Scholar 
59.Ackerly, D. & Cornwell, W. A trait-based approach to community assembly: Partitioning of species trait values into within- and among-community components. Ecol. Lett. 10, 135–145 (2007).CAS 
PubMed 
Article 

Google Scholar 
60.Rijkers, T., Pons, T. L. & Bongers, F. The effect of tree height and light availability on photosynthetic leaf traits of four neotropical species differing in shade tolerance. Funct. Ecol. 14, 77–86 (2000).Article 

Google Scholar 
61.de Bello, F. et al. Partitioning of functional diversity reveals the scale and extent of trait convergence and divergence. J. Veg. Sci. 20, 475–486 (2009).Article 

Google Scholar 
62.Ding, Y., Zang, R., Lu, X. & Huang, J. The impacts of selective logging and clear-cutting on woody plant diversity after 40years of natural recovery in a tropical montane rain forest, south China. Sci. Total Environ. 579, 1683–1691 (2017).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar  More

We now derive the metapopulation model used in this paper. We start by deriving a general metapopulation model that is based on the seminal work of Levin50. Assuming that the inter-patch migrations are detailed-balanced, we make use of the formulation in Eq. (8) to derive a balanced metapopulation model. We then show that the balanced model admits a unique coexistence equilibrium that is asymptotically stable if the dispersal network is heterogeneous, whereas the same equilibrium is neutrally stable in the case of a homogeneous network.General metapopulation modelMathematical models based on traditional metapopulation theory usually assume that the metapopulation is made up of many neighboring spatially homogeneous habitat patches connected via dispersal. Consider an interconnected network of m discrete patches each being inhabited by the same n species. In addition, assume that species can migrate from one patch to some or all of the other patches. The rate of migration of each species between two patches is directly proportional to the proportion of the particular species in the originating patch, with a (nonnegative) constant of proportionality being the same across species. This constant of proportionality will be referred to as the rate constant associated with the migration. It is assumed that if there is migration between two given patches, then it is bidirectional, i.e., the rate constant of migration from j to k is strictly positive if and only the same holds for the migration from k to j. Just like in the case of a reversible single-species chemical reaction network, inter-patch migrations may be described by a weighted symmetric directed graph (G_2=(V_2,E_2)) where (V_2={1,ldots ,m}) is the set of patches (vertices) and an edge ((j,k)in E_2) means that every species can migrate from patch j to patch k. Finally, it is also assumed that the graph (G_2) corresponding to the inter-patch migration is connected, i.e., there is a path between every two distinct vertices of the graph.The flow of species between the patches can be summarized in a weighted (mtimes m) adjacency matrix ({mathbf {A}}) with entry (A_{jk}) being equal to the rate constant of migration of species from the (j{ {text {th}}}) to the (k{ {text {th}}}) patch. The diagonal elements of ({mathbf {A}}) are hence equal to 0. Due to the bidirectional nature of migration, it holds that (A_{jk} >0 Leftrightarrow A_{kj} >0) and (A_{jk}=0 Leftrightarrow A_{kj}=0), for any (jne k). Let (Delta =text {diag}(delta _1,ldots ,delta _m)) denote the m-dimensional diagonal matrix whose (j{ {text {th}}}) entry is given by$$begin{aligned} delta _{j}=sum _{k=1}^{m} A_{jk}. end{aligned}$$Define ({mathbf {L}}:=Delta -{mathbf {A}}^top ). Note that$$begin{aligned} (mathbb {1}^m)^{top }{mathbf {L}}=(mathbb {1}^m)^{top }Delta -big ({mathbf {A}}mathbb {1}^mbig )^{top }=({mathbf {0}}^m)^top . end{aligned}$$Let ({mathbf {x}}in S^{mn}), with (x_{i,j}) the proportion of species i in patch j across the entire metapopulation, then the net migration rate (psi _{i,j}) of species i from other patches to patch j is given by$$begin{aligned} psi _{i,j}=sum _{k=1}^{m}A_{kj}x_{i,k}-sum _{k=1}^{m}A_{jk}x_{i,j}=sum _{k=1}^{m}A_{kj}x_{i,k}-delta _{j}x_{i,j}=-sum _{k=1}^{m}L_{jk}x_{i,k}. end{aligned}$$Let us denote (Psi _i:=left( psi _{i,1},psi _{i,2},ldots ,psi _{i,m}right) ^{top }) and ({mathbf {r}}_i:=left( x_{i,1},x_{i,2},ldots ,x_{i,m}right) ^{top }), then$$begin{aligned} Psi _{i}=-{mathbf {L}}{mathbf {r}}_{i}. end{aligned}$$
(9)
Within each patch, the proportions of species are affected by other patches only via migration. Let (phi _{i,j}) denote the rate of change of the proportion of species i in patch j in the absence of migration. Since the dominance relationships among the species (described by a tournament matrix ({mathbf {T}})) are assumed to be the same for all patches and since the habitat patches are spatially homogeneous, the expression for (phi _{i,j}) is given by the right-hand side of System (1):$$begin{aligned} phi _{i,j}=x_{i,j}left( {mathbf {T}}{mathbf {p}}_{j}right) _{i}, end{aligned}$$
(10)
where ({mathbf {p}}_j:=left( x_{1,j},x_{2,j},ldots , x_{n,j}right) ^{top }), (i=1,ldots ,n) and (j=1,ldots ,m). Assuming migration among the patches, the proportion of a species within a patch is influenced by two factors: the first is the interaction with other species within the patch and the second is the migration of that particular species to or from other patches. Thus, the metapopulation model describing the dynamics of the n species in the m-patch network is described by the system of mn differential equations;$$begin{aligned} {dot{x}}_{i,j}=phi _{i,j}+psi _{i,j}=x_{i,j}left( {mathbf {T}}{mathbf {p}}_{j}right) _{i}-left( {mathbf {L}}{mathbf {r}}_{i}right) _{j},qquad i=1,ldots ,n,quad j=1,ldots ,m . end{aligned}$$
(11)
This system evolves on the unit simplex (S^{mn}).
Proposition 2

The unit simplex (S^{mn}) is positively invariant for System (11).

Proof
To show the invariance of the unit simplex (S^{mn}) under the flow of System (11), it suffices to show that each of the faces of the simplex cannot be crossed, i.e., the vector field points inward from the faces of (S^{mn}).
On the one hand, if (x_{i,j}=0) for some i, j, then$$begin{aligned}{dot{x}}_{i,j}=sum _{k=1}^{m}A_{kj}x_{i,k}ge 0,end{aligned}$$which implies that (x_{i,j}=0) cannot be crossed from positive to negative. In an ecological context, this condition simply states the obvious fact that an extinct species is in no danger of declining. On the other hand, if (x_{i,j}=1) for some i, j, then obviously (x_{l,k}=0) for any (lne i) or (kne j) and$$begin{aligned}{dot{x}}_{i,j}=-delta _{j}< 0.end{aligned}$$Hence, the vector field associated with System (11) points inward from the faces of (S^{mn}). So, (S^{mn}) is positively invariant under the flow of System (11). (square ) Note that Proposition 2 does not exclude the solution trajectories of System (11) from approaching the boundary equilibria of the system as (trightarrow infty ). We call metapopulation model (11) persistent if for every ({mathbf {x}}_0in S^{mn}_{+}), the (omega )-limit set (omega ({mathbf {x}}_0)) does not intersect the boundary of (S^{mn}). In other words, a metapopulation model is persistent if the initial existence of all the species implies that none of the species goes extinct with the passage of time.Balanced homogeneous and heterogeneous metapopulation modelsWe say that the inter-patch migration of a metapopulation model is detailed balanced if the overall migration rate of any species between any two patches is zero for a certain positive set of proportions of that species in the different patches. From the theory of detailed-balanced reaction networks described in “Detailed-balanced single species mass action reaction networks” section, it follows that a detailed-balanced inter-patch migration network corresponds to a detailed-balanced single species mass action reaction network. Let B denote the incidence matrix corresponding to the directed graph (G_2) describing the inter-patch migrations and let r denote the number of edges in (G_2). Comparing Eqs. (8) and (9), it follows that if the inter-patch migration is detailed balanced, then there exist diagonal matrices ({mathcal {K}}in {mathbb {R}}^{rtimes r}) and ({mathbf {Z}}^*in {mathbb {R}}^{mtimes m}) with positive diagonal entries such that ((mathbb {1}^m)^{top }{mathbf {Z}}^*mathbb {1}^m=1) and$$begin{aligned} {mathbf {L}}={mathbf {B}}{mathcal {K}}{mathbf {B}}^{top }({mathbf {Z}}^*)^{-1}. end{aligned}$$Let ({mathbf {Z}}^*=text { diag}({mathbf {z}}^*)). Equation (9) can now be rewritten as$$begin{aligned} Psi _{i}=-{mathbf {B}}{mathcal {K}}{mathbf {B}}^{top }left( frac{{mathbf {r}}_{i}}{{mathbf {z}}^{*}}right) . end{aligned}$$ (12) Henceforth in this manuscript, we restrict our analysis to metapopulation models of type (11) for which the interactions within each patch correspond to a tournament with a completely mixed optimal strategy and whose inter-patch migration is detailed balanced. Such metapopulation models will be referred to as balanced metapopulation models.We have seen earlier in “Species interactions and tournament matrices” section that if the interactions within every patch correspond to a tournament with a completely mixed optimal strategy, then the corresponding mean-field model admits a unique coexistence equilibrium ({mathbf {y}}^*in S^{n}_{+}) with ({mathbf {T}}{mathbf {y}}^*={mathbf {0}}^n). Thus, for a balanced metapopulation model, System (10) can be rewritten as$$begin{aligned} phi _{i,j}=x_{i,j}left( mathbf {TY}^*left( frac{{mathbf {p}}_{j}}{{mathbf {y}}^*}right) right) _i, end{aligned}$$ (13) where ({mathbf {Y}}^*:=) diag(({mathbf {y}}^*)). Consequently, from Eqs. (11)–(13), it follows that the dynamics of a balanced metapopulation model containing n species and m patches can be described by mn differential equations$$begin{aligned} {dot{x}}_{i,j}=x_{i,j}left( mathbf {TY}^*left( frac{{mathbf {p}}_{j}}{{mathbf {y}}^*}right) right) _i-left( {mathbf {B}}{mathcal {K}}{mathbf {B}}^{top }left( frac{{mathbf {r}}_{i}}{{mathbf {z}}^{*}}right) right) _{j} ,qquad i=1,ldots ,n,quad j=1,ldots ,m . end{aligned}$$ (14) If all the elements of ({mathbf {z}}^*) in the above equation are equal, i.e., if (z_j^*=frac{1}{m}) for (j=1,ldots ,m), then we say that the balanced metapopulation model is homogeneous, otherwise we call it heterogeneous. Whether a balanced metapopulation model is homogeneous or not can be checked from the adjacency matrix ({mathbf {A}}) corresponding to its inter-patch migration graph (G_2). If ({mathbf {A}}) is symmetric, then the model is homogeneous, otherwise it is heterogeneous. Remark 3 In35, the authors assume that migrations from one patch to other patches are random with a probability of migration (or migration constant) equal to the reciprocal of the number of dispersal links from a patch to other patches. They thus define a dispersal graph to be homogeneous if all nodes have the same degree (number of links), otherwise the graph is heterogeneous. With this definition, homogeneity, in general, is equivalent to the existence of cycles in the dispersal graph, whereas heterogeneity is equivalent to their absence. However, with our new definition, it is clear that this is not necessary. An example of such a case is shown in Fig. 2.Figure 2Left: A heterogeneous dispersal graph according to35. Right: A homogeneous dispersal graph according to our definition.Full size image Coexistence equilibrium and its uniquenessIn this section, we present a theorem that gives an expression for a coexistence equilibrium of a balanced metapopulation model. Before we state our main theorem in this section, we need the following lemma. Lemma 4 Let ({mathbf {B}}in {mathbb {R}}^{mtimes r}) denote the incidence matrix of a finite connected directed graph (G_2) and let ({mathcal {K}}in {mathbb {R}}^{rtimes r}) denote a diagonal matrix with positive diagonal entries. For any ({mathbf {w}}in {mathbb {R}}_+^m), it holds that (-{mathbf {w}}^{top }{mathbf {B}}{mathcal {K}}{mathbf {B}}^{top }left( frac{mathbb {1}^m}{{mathbf {w}}}right) ge 0). Moreover (-{mathbf {w}}^{top }{mathbf {B}}{mathcal {K}}{mathbf {B}}^{top }left( frac{mathbb {1}^m}{{mathbf {w}}}right) = 0) if and only if ({mathbf {w}}=qmathbb {1}^m), where (qin {mathbb {R}}_+). Proof Assume that the (p{ {text {th}}}) edge of the graph (G_2) is directed from vertex (i_p) to vertex (j_p). Hence, (B_{i_pp}=-1), (B_{j_pp}=1) and (B_{kp}=0) for (i_pne kne j_p). Thus,$$begin{aligned} -{mathbf {w}}^{top }{mathbf {B}}{mathcal {K}}{mathbf {B}}^{top }left( frac{mathbb {1}^m}{{mathbf {w}}}right) =sum _{p=1}^m(w_{j_p}-w_{i_p})kappa _pleft( frac{1}{w_{i_p}}-frac{1}{w_{j_p}}right) =sum _{p=1}^mfrac{kappa _p}{w_{i_p}w_{j_p}}left( w_{j_p}-w_{i_p}right) ^2ge 0. end{aligned}$$Moreover, (-{mathbf {w}}^{top }{mathbf {B}}{mathcal {K}}{mathbf {B}}^{top }left( frac{mathbb {1}^m}{{mathbf {w}}}right) =0) if and only if (w_{j_p}=w_{i_p}) for (p=1,ldots ,m), which is equivalent with ({mathbf {B}}^{top }{mathbf {w}}={mathbf {0}}^r). Since the graph (G_2) is connected, we recall from48 that (text {rank}({mathbf {B}})=m-1) and furthermore (text {ker}({mathbf {B}}^{top })=mathbb {1}^m). Therefore ({mathbf {B}}^{top }{mathbf {w}}={mathbf {0}}^r) if and only if ({mathbf {w}}=qmathbb {1}^m), where (qin {mathbb {R}}_+). This completes the proof. (square ) We now state the main theorem of this section. Theorem 5 A balanced metapopulation model described by System (14) admits a unique coexistence equilibrium ({mathbf {x}}^*in S^{mn}_{+}). The proportion (x_{i,j}^{*}) of species i in patch j at the unique coexistence equilibrium is given by$$begin{aligned} x_{i,j}^*=y^{*}_iz^{*}_j. end{aligned}$$ (15) for (i=1,ldots ,n) and (j=1,ldots ,m). Proof We divide the proof into two parts. In the first part we prove that System (15) indeed yields an equilibrium for the model. In the second part, we prove that this coexistence equilibrium is unique. Let us define$$begin{aligned} {mathbf {p}}_{j}^*:=left( x_{1,j}^*, x_{2,j}^*, ldots , x_{n,j}^*right) ^top =z_j^*{mathbf {y}}^*; quad {mathbf {r}}_{i}^{*}:=left( x_{i,1}^{*}, x_{i,2}^{*}, ldots , x_{i,m}^{*}right) ^top =y_i^*{mathbf {z}}^*. end{aligned}$$For ({mathbf {x}}^*) to be an equilibrium of System (14), it should render the right-hand side equal to zero. Note that$$begin{aligned} mathbf {TY}^*left( frac{{mathbf {p}}_{j}^*}{{mathbf {y}}^*}right) =z_j^*mathbf {TY}^*mathbb {1}^n=z_j^*{mathbf {T}}{mathbf {y}}^{*}={mathbf {0}}^n end{aligned}$$and$$begin{aligned} {mathbf {B}}{mathcal {K}}{mathbf {B}}^{top }left( frac{{mathbf {r}}_{i}^*}{{mathbf {z}}^{*}}right) =y_i^*{mathbf {B}}{mathcal {K}}{mathbf {B}}^{top }mathbb {1}^m={mathbf {0}}^m. end{aligned}$$In addition,$$begin{aligned} (mathbb {1}^{mn})^{top }{mathbf {x}}^*=sum _{i=1}^nsum _{j=1}^mx_{i,j}^*=sum _{i=1}^{n}y_i^*sum _{j=1}^mz_j^{*}=1. end{aligned}$$Thus, ({mathbf {x}}^*) is a coexistence equilibrium of System (14). Assume that there exists another coexistence equilibrium ({mathbf {x}}^{**}in , S^{mn}_{+}). Let (x_{i,j}^{**}) denote the corresponding proportion of species i in patch j and define$$begin{aligned} {mathbf {p}}_{j}^{**}:=left( x_{1,j}^{**}, x_{2,j}^{**}, ldots , x_{n,j}^{**}right) ^top ; qquad {mathbf {r}}_{i}^{**}:=left( x_{i,1}^{**}, x_{i,2}^{**}, ldots , x_{i,m}^{**}right) ^top . end{aligned}$$It follows that for any i, j it holds that$$begin{aligned} x_{i,j}^{**}left( mathbf {TY}^*left( frac{{mathbf {p}}_{j}^{**}}{{mathbf {y}}^*}right) right) _i-left( {mathbf {B}}{mathcal {K}}{mathbf {B}}^{top }left( frac{{mathbf {r}}_{i}^{**}}{{mathbf {z}}^{*}}right) right) _{j}=0. end{aligned}$$ (16) Multiplying both sides of this equality with (frac{x_{i,j}^*}{x_{i,j}^{**}}), we get$$begin{aligned} x_{i,j}^*left( mathbf {TY}^*left( frac{{mathbf {p}}_{j}^{**}}{{mathbf {y}}^*}right) right) _i- frac{x_{i,j}^*}{x_{i,j}^{**}} left( {mathbf {B}}{mathcal {K}}{mathbf {B}}^{top }left( frac{{mathbf {r}}_{i}^{**}}{{mathbf {z}}^{*}}right) right) _{j}=0. end{aligned}$$Summing the left-hand side of the above expression over the different species and patches, we get$$begin{aligned} sum _{j=1}^msum _{i=1}^nx_{i,j}^*left( mathbf {TY}^*left( frac{{mathbf {p}}_{j}^{**}}{{mathbf {y}}^*}right) right) _i- sum _{i=1}^nsum _{j=1}^mfrac{x_{i,j}^*}{x_{i,j}^{**}} left( {mathbf {B}}{mathcal {K}}{mathbf {B}}^{top }left( frac{{mathbf {r}}_{i}^{**}}{{mathbf {z}}^{*}}right) right) _{j}=0. end{aligned}$$ (17) Now consider the two terms in the left-hand side of the above equality separately. For the first term, note that for any j it holds that$$begin{aligned} sum _{i=1}^nx_{i,j}^*left( mathbf {TY}^*left( frac{{mathbf {p}}_{j}^{**}}{{mathbf {y}}^*}right) right) _i= & {} sum _{i=1}^nx_{i,j}^*left( {mathbf {T}}{mathbf {p}}_{j}^{**}right) _{i} = sum _{i=1}^{n}x_{i,j}^{*}left( sum _{l=1}^{n}T_{il}x_{l,j}^{**}right) =-sum _{l=1}^{n}x_{l,j}^{**}left( sum _{i=1}^{n}T_{li}x_{i,j}^{*}right) \= & {} -sum _{l=1}^{n}x_{l,j}^{**}left( sum _{i=1}^nT_{li}y_i^{*}z_j^*right) =-z_j^*sum _{l=1}^{n}x_{l,j}^{**}({mathbf {T}}{mathbf {y}}^*)_l=0. end{aligned}$$Hence,$$begin{aligned} sum _{j=1}^msum _{i=1}^nx_{i,j}^*left( mathbf {TY}^*left( frac{{mathbf {p}}_{j}^{**}}{{mathbf {y}}^*}right) right) _i=0. end{aligned}$$For the second term, we find$$begin{aligned} -sum _{i=1}^nsum _{j=1}^mfrac{x_{i,j}^*}{x_{i,j}^{**}} left( {mathbf {B}}{mathcal {K}}{mathbf {B}}^{top }left( frac{{mathbf {r}}_{i}^{**}}{{mathbf {z}}^{*}}right) right) _{j}=-sum _{i=1}^ny_i^*sum _{j=1}^mfrac{z_j^*}{x_{i,j}^{**}}left( {mathbf {B}}{mathcal {K}}{mathbf {B}}^{top }left( frac{{mathbf {r}}_{i}^{**}}{{mathbf {z}}^{*}}right) right) _{j}=-sum _{i=1}^ny_i^*left( frac{{mathbf {z}}^{*}}{{mathbf {r}}_{i}^{**}}right) ^{top }{mathbf {B}}{mathcal {K}}{mathbf {B}}^{top }left( frac{{mathbf {r}}_{i}^{**}}{{mathbf {z}}^{*}}right) . end{aligned}$$Thus, Eq. (17) can be simplified as$$begin{aligned} -sum _{i=1}^ny_i^*left( frac{{mathbf {z}}^{*}}{{mathbf {r}}_{i}^{**}}right) ^{top }{mathbf {B}}{mathcal {K}}{mathbf {B}}^{top }left( frac{{mathbf {r}}_{i}^{**}}{{mathbf {z}}^{*}}right) =0. end{aligned}$$Since (y_i^* >0) for (i=1,ldots ,n), it holds for any (i=1,ldots ,n) that$$begin{aligned} -left( frac{{mathbf {z}}^{*}}{{mathbf {r}}_{i}^{**}}right) ^{top }{mathbf {B}}{mathcal {K}}{mathbf {B}}^{top }left( frac{{mathbf {r}}_{i}^{**}}{{mathbf {z}}^{*}}right) =0. end{aligned}$$
(18)
From Eq. (18) and Lemma 4, it follows that ({mathbf {r}}_{i}^{**}=q_i{mathbf {z}}^*) with (q_iin {mathbb {R}}_+) for (i=1,ldots ,n). Thus, (x_{i,j}^{**}=q_iz_{j}^*) and ({mathbf {p}}_{j}^{**}=z_j^*{mathbf {q}}) for (i=1,ldots ,n) and (j=1,ldots ,m). Substituting the latter in the left-hand side of Eq. (16), we get$$begin{aligned} x_{i,j}^{**}left( mathbf {TY}^*left( frac{{mathbf {p}}_{j}^{**}}{{mathbf {y}}^*}right) right) _i-left( {mathbf {B}}{mathcal {K}}{mathbf {B}}^{top }left( frac{{mathbf {r}}_{i}^{**}}{{mathbf {z}}^{*}}right) right) _{j}=q_i{z_j^*}^2left( {mathbf {T}}{mathbf {Y}}^*left( frac{{mathbf {q}}}{{mathbf {y}}^*}right) right) _i-q_ileft( {mathbf {B}}{mathcal {K}}{mathbf {B}}^{top }mathbb {1}^mright) _j=q_i{z_j^*}^2(mathbf {Tq})_i. end{aligned}$$Since (q_i >0) for (i=1,ldots ,n), for Eq. (16) to hold, we should have (mathbf {Tq}={mathbf {0}}^n). Also note that$$begin{aligned} (mathbb {1}^{mn})^{top }{mathbf {x}}^{**}=sum _{i=1}^nsum _{j=1}^{m}x_{i,j}^{**}=sum _{i=1}^nq_isum _{j=1}^mz_j^* =sum _{i=1}^nq_i=1. end{aligned}$$Since the metapopulation model is balanced, it follows that ({mathbf {q}}={mathbf {y}}^*). Thus, (x_{i,j}^{**}=y_i^*z_j^*=x_{i,j}^*) for (i=1,ldots ,n) and (j=1,ldots ,m). This proves the uniqueness of the coexistence equilibrium ({mathbf {x}}^*). (square )
We now give examples of two balanced metapopulation models.

Example 1

It is easy to verify that the network shown in Fig. 3 corresponds to a balanced metapopulation model governed by System (14) with$$begin{aligned} {mathbf {T}} = left[ begin{array}{rrr} 0 &{}quad 1 &{}quad -1\ -1 &{}quad 0 &{}quad 1\ 1 &{}quad -1 &{}quad 0 end{array}right] ; quad {mathbf {B}} = left[ begin{array}{rrr} -1 &{}quad 0 &{}quad 1\ 1 &{}quad -1 &{}quad 0\ 0 &{}quad 1 &{}quad -1 end{array}right] ; end{aligned}$$({mathbf {y}}^*=left( frac{1}{3}, frac{1}{3}, frac{1}{3} right) ^{top }), ({mathbf {z}}^*=left( frac{1}{5}, frac{2}{5}, frac{2}{5} right) ^{top }) and ({mathcal {K}}=text { diag}left( frac{1}{10},frac{3}{10},frac{1}{10}right) ). Note that this metapopulation model is heterogeneous. From Theorem 5, it follows that the species proportions at the unique coexistence equilibrium for this model are given by (x_{i,1}^*=frac{1}{15}) and (x_{i,2}^*=x_{i,3}^*= frac{2}{15}).Figure 3A metapopulation network composed of three patches. Each patch contains a local population composed of three species (1, 2 and 3), in cyclic competition, as shown by the black arrows. The red arrows denote migrations among the patches in the directions shown.Full size image

Example 2

It is easy to verify that the network shown in Fig. 4 corresponds to a balanced metapopulation model governed by System (14) with$$begin{aligned} {mathbf {T}} = left[ begin{array}{rrr} 0 &{}quad 1 &{} quad -1\ -1 &{}quad 0 &{} quad 1\ 1 &{}quad -1 &{}quad 0 end{array}right] ; quad {mathbf {B}} = left[ begin{array}{rrr} 1 &{}quad -1\ 0 &{}quad 1\ -1 &{}quad 0 end{array}right] ; end{aligned}$$({mathbf {y}}^{*}={mathbf {z}}^*=left( frac{1}{3}, frac{1}{3}, frac{1}{3} right) ^{top }) and ({mathcal {K}}=frac{1}{3}text { diag}(mathbb {1}_2)). Note that this metapopulation model is homogeneous. From Theorem 5, it follows that the species proportions at the unique coexistence equilibrium in this case are all given by (x_{i,j}^*=frac{1}{9}) for (i,j= 1,2,3).Figure 4A metapopulation network composed of three patches. Species can migrate from patch 1 to the other two patches and vice versa. However, there exists no migrations between patches 2 and 3.Full size image
StabilityWe now prove the local stability of the unique coexistence equilibrium corresponding to the balanced metapopulation model (14). For the proof, we make use of the same Lyapunov function as in “Neutral stability” section, coupled with LaSalle’s invariance principle51, (52, Section 4.2), (53, pp. 188–189).

Theorem 6

Consider the balanced metapopulation model (14) with coexistence equilibrium ({mathbf {x}}^*in , S^{mn}_{+}).

1.

If the model is heterogeneous, then ({mathbf {x}}^*) is locally asymptotically stable w.r.t. all initial conditions in (S^{mn}_{+}) in the neighbourhood of ({mathbf {x}}^*). Furthermore, if the model is persistent, then ({mathbf {x}}^*) is globally asymptotically stable w.r.t. all initial conditions in (S^{mn}_{+}).

2.

If the model is homogeneous and persistent, then as (trightarrow infty ), the solution trajectories converge to a limit cycle satisfying the equation ({dot{x}}_{i,j}=x_{i,j}({mathbf {T}}{mathbf {p}}_{j})_i) with (x_{i,j}=x_{i,k}), for (i=1,ldots ,n) and (j,k=1,ldots ,m).

Proof
Let (x_{i,j}) denote the proportion of species i in patch j. Assuming that ({mathbf {x}}in S^{mn}_{+}), consider the Lyapunov function$$begin{aligned} V({mathbf {x}})=-(mathbf {x^{*}})^{top }text {Ln}left( frac{{mathbf {x}}}{{mathbf {x}}^*}right) . end{aligned}$$
(19)
By Gibbs inequality, V(x) is positive on (S^{mn}_{+}) and is equal to zero only if ({mathbf {x}}={mathbf {x}}^*). Taking the time derivative of V, we have$$begin{aligned} {dot{V}}({mathbf {x}})=-sum _{j=1}^msum _{i=1}^{n}left( frac{x_{i,j}^{*}}{x_{i,j}}right) {dot{x}}_{i,j}. end{aligned}$$From Eq. (14), it follows that$$begin{aligned} {dot{V}}({mathbf {x}})= -sum _{j=1}^msum _{i=1}^nx_{i,j}^*left( mathbf {TY}^*left( frac{{mathbf {p}}_{j}}{{mathbf {y}}^*}right) right) _i+ sum _{i=1}^nsum _{j=1}^mfrac{x_{i,j}^*}{x_{i,j}} left( {mathbf {B}}{mathcal {K}}{mathbf {B}}^{top }left( frac{{mathbf {r}}_{i}}{{mathbf {z}}^{*}}right) right) _{j}. end{aligned}$$As in the proof of Theorem 5, it can be verified that$$begin{aligned} sum _{j=1}^msum _{i=1}^nx_{i,j}^*left( mathbf {TY}^*left( frac{{mathbf {p}}_{j}}{{mathbf {y}}^*}right) right) _i=0 end{aligned}$$and$$begin{aligned} sum _{i=1}^nsum _{j=1}^mfrac{x_{i,j}^*}{x_{i,j}} left( {mathbf {B}}{mathcal {K}}{mathbf {B}}^{top }left( frac{{mathbf {r}}_{i}}{{mathbf {z}}^{*}}right) right) _{j}=sum _{i=1}^ny_i^*left( frac{{mathbf {z}}^{*}}{{mathbf {r}}_{i}}right) ^{top }{mathbf {B}}{mathcal {K}}{mathbf {B}}^{top }left( frac{{mathbf {r}}_{i}}{{mathbf {z}}^{*}}right) . end{aligned}$$Thus,$$begin{aligned} {dot{V}}({mathbf {x}})=sum _{i=1}^ny_i^*left( frac{{mathbf {z}}^{*}}{{mathbf {r}}_{i}}right) ^{top }{mathbf {B}}{mathcal {K}}{mathbf {B}}^{top }left( frac{{mathbf {r}}_{i}}{{mathbf {z}}^{*}}right) . end{aligned}$$Since (y_i^* >0) for (i=1,ldots ,n), it follows from Lemma 4 that ({dot{V}}({mathbf {x}})le 0) and ({dot{V}}({mathbf {x}})=0) if and only if ({mathbf {r}}_i=q_i{mathbf {z}}^*) with (q_iin {mathbb {R}}_+), for (i=1,ldots ,n). Thus,$$begin{aligned} x_{i,j}=q_iz_j^*, end{aligned}$$
(20)
for (i= 1,ldots ,n) and (j=1,ldots ,m). Since ((mathbb {1}^{mn})^{top }{mathbf {x}}=1), we obtain$$begin{aligned} sum _{i=1}^nsum _{j=1}^{m}x_{i,j}=sum _{i=1}^nq_isum _{j=1}^mz_j^* =sum _{i=1}^nq_i=1. end{aligned}$$Let ({mathcal {E}}subset S^{mn}_{+}) be the set of all vectors ({mathbf {x}}) for which condition (20) is satisfied with ((mathbb {1}^n)^{top }{mathbf {q}}=1). We now determine the largest subset of ({mathcal {E}}) that is positively invariant w.r.t. System (14). Assume that ({mathbf {x}}) continuously takes values from ({mathcal {E}}) and satisfies System (14). Since ({mathbf {x}}) takes values from ({mathcal {E}}), we have ({dot{x}}_{i,j}=z_j^*{dot{q}}_i). Since ({mathbf {x}}) also satisfies System (14), we have$$begin{aligned} {dot{x}}_{i,j}=x_{i,j}left( {mathbf {T}}{mathbf {p}}_{j}right) _i-left( {mathbf {B}}{mathcal {K}}{mathbf {B}}^{top }left( frac{{mathbf {r}}_{i}^{*}}{{mathbf {z}}^{*}}right) right) _{j}=q_i{z_j^*}^2(mathbf {Tq})_i-q_ileft( {mathbf {B}}{mathcal {K}}{mathbf {B}}^{top }mathbb {1}^mright) _j=q_i{z_j^*}^2(mathbf {Tq})_i. end{aligned}$$Thus, (z_j^*{dot{q}}_i=q_i{z_j^*}^2(mathbf {Tq})_i) which implies that$$begin{aligned} {dot{q}}_i=z_j^*q_i(mathbf {Tq})_i, end{aligned}$$
(21)
for (i=1,ldots ,n) and (j=1,ldots ,m). We now consider two cases.
Case 1: The model is heterogeneous, i.e., the vector ({mathbf {z}}^*) is not parallel to (mathbb {1}^m).
In this case, Eq. (21) will be satisfied only if (q_i(mathbf {Tq})_i=0) for (i=1,ldots ,n). Since (q_iin {mathbb {R}}_+) for (i=1,ldots ,n), it follows that (mathbf {Tq}={mathbf {0}}^n). Since ((mathbb {1}^n)^{top }{mathbf {q}}=1), we have ({mathbf {q}}={mathbf {y}}^*). This implies that (x_{i,j}=y_i^*z_j^*=x_{i,j}^*) for (i=1,ldots ,n) and (j= 1,ldots ,m). Thus, the largest subset of ({mathcal {E}}) that is positively invariant w.r.t. System (14) consists of just the unique equilibrium ({mathbf {x}}^*in S^{mn}_{+}). By LaSalle’s invariance principle, it follows that the equilibrium ({mathbf {x}}^*) is locally asymptotically stable w.r.t. all initial conditions in (S^{mn}_{+}) in the neighbourhood of ({mathbf {x}}^*), and globally asymptotically stable w.r.t. all initial conditions in (S^{mn}_{+}) provided that System (14) is persistent.

Case 2: The model is homogeneous, i.e.
({mathbf {z}}^*=frac{1}{m}mathbb {1}^m)

In this case, Eq. (21) takes the form ({dot{q}}_i=frac{q_i}{m}(mathbf {Tq})_i). We have (x_{i,j}=q_iz_j^*=frac{q_i}{m}) and$$begin{aligned} {dot{x}}_{i,j}=frac{{dot{q}}_i}{m}=frac{q_i}{m^2}(mathbf {Tq})_i=x_{i,j}({mathbf {T}}{mathbf {p}}_{j})_i. end{aligned}$$Consequently, the largest subset of ({mathcal {E}}) that is positively invariant w.r.t. System (14) consists of all vectors ({mathbf {x}}(t)in , S^{mn}_{+}) satisfying ({dot{x}}_{i,j}=x_{i,j}({mathbf {T}}{mathbf {p}}_{j})_i) with (x_{i,j}=x_{i,k}) for (i=1,ldots ,n) and (j,k=1,ldots ,m). The proof for Case 2 again follows from LaSalle’s invariance principle. (square )
The above results can be illustrated by simulating System (14) for the metapopulation models shown in Fig. 3 and 4 in Examples 1 and 2, respectively. The results of the simulations are shown in Figs. 5 and 6, respectively.Figure 5Left: Dynamics of the metapopulation model in Fig. 3 for patches 1 and 3 showing asymptotic stability of the coexistence equilibrium. Right: The time evolution of the proportion of species 1 in the three patches.Full size imageFigure 6Left: Dynamics of the metapopulation model in Fig. 4 for patches 1 and 3 showing a limit cycle arising from the neutral stability of the coexistence equilibrium. Right: Time evolution of the proportion of species 1 in the three patches. Note that the dynamics in all patches are the same and thus the three graphs overlap.Full size image More

1.Jackson, C. & Robertson, M. Predicting the potential distribution of an endangered cryptic subterranean mammal from few occurrence records. J. Nat. Conserv. https://doi.org/10.1016/j.jnc.2010.06.006 (2011).Article 

Google Scholar 
2.Rondinini, C. et al. Global habitat suitability models of terrestrial mammals. Philos. Trans. R. Soc. B Biol. Sci. 366, 2633–2641 (2011).Article 

Google Scholar 
3.Guisan, A. & Thuiller, W. Predicting species distribution: Offering more than simple habitat models. Ecol. Lett. 8, 993–1009 (2005).Article 

Google Scholar 
4.Yang, X.-Q., Kushwaha, S. P. S., Saran, S., Xu, J. & Roy, P. S. Maxent modeling for predicting the potential distribution of medicinal plant, Justicia adhatoda L. in Lesser Himalayan foothills. Ecol. Eng. 51, 83–87 (2013).CAS 
Article 

Google Scholar 
5.Ouyang, Z., Liu, J., Xiao, H., Tan, Y. & Zhang, H. An assessment of giant panda habitat in Wolong Nature Reserve. Acta Ecol. Sin. 11, 1869–1874 (2001).
Google Scholar 
6.Schadt, S. et al. Assessing the suitability of central European landscapes for the reintroduction of Eurasian lynx. J. Appl. Ecol. 39, 189–203 (2002).Article 

Google Scholar 
7.Su, J., Aryal, A., Nan, Z. & Ji, W. Climate change-induced range expansion of a subterranean rodent: Implications for rangeland management in Qinghai-Tibetan Plateau. PLoS One 10, e0138969 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
8.Srivastava, V., Griess, V. C. & Padalia, H. Mapping invasion potential using ensemble modelling. A case study on Yushania maling in the Darjeeling Himalayas. Ecol. Model. 385, 35–44 (2018).Article 

Google Scholar 
9.Phillips, S. J., Anderson, R. P. & Schapire, R. E. Maximum entropy modeling of species geographic distributions. Ecol. Model. 190, 231–259 (2006).Article 

Google Scholar 
10.Raffini, F. et al. From nucleotides to satellite imagery: Approaches to identify and manage the invasive pathogen Xylella fastidiosa and its insect vectors in Europe. Sustainability 12, 4508 (2020).CAS 
Article 

Google Scholar 
11.Clements, G. R. et al. Predicting the distribution of the Asian Tapir (Tapirus indicus) in Peninsular Malaysia using maximum entropy modelling. Integr. Zool. 7, 400–406 (2012).PubMed 
Article 

Google Scholar 
12.Hijmans, R. J. & Graham, C. H. The ability of climate envelope models to predict the effect of climate change on species distributions. Glob. Chang. Biol. 12, 2272–2281 (2006).ADS 
Article 

Google Scholar 
13.Phillips, S. J. & Dudík, M. Modeling of species distributions with Maxent: New extensions and a comprehensive evaluation. Ecography (Cop.) 31, 161–175 (2008).Article 

Google Scholar 
14.Cassini, M. H. Ecological principles of species distribution models: The habitat matching rule. 2057–2065. https://doi.org/10.1111/j.1365-2699.2011.02552.x (2011).15.Elith, J. & Leathwick, J. R. Species distribution models: Ecological explanation and prediction across space and time. Annu. Rev. Ecol. Evol. Syst. 40, 677–697 (2009).Article 

Google Scholar 
16.Mac Nally, R. Regression and model-building in conservation biology, biogeography and ecology: The distinction between–and reconciliation of–‘predictive’ and ‘explanatory’models. Biodivers. Conserv. 9, 655–671 (2000).Article 

Google Scholar 
17.Jaynes, E. T. Information theory and statistical mechanics. II. Phys. Rev. 108, 171–190 (1957).ADS 
MathSciNet 
MATH 
Article 

Google Scholar 
18.Jaynes, E. T. Probability Theory as Logic BT – Maximum Entropy and Bayesian Methods. In (ed. Fougère, P. F.) 1–16 (Springer, Netherlands, 1990). https://doi.org/10.1007/978-94-009-0683-9_1.19.Jaynes, E. T. Probability Theory: The Logic of Science (Cambridge University Press, Cambridge, 2003).MATH 
Book 

Google Scholar 
20.Stockwell, D. The GARP modelling system: Problems and solutions to automated spatial prediction. Int. J. Geogr. Inf. Sci. 13, 143–158 (1999).Article 

Google Scholar 
21.Townsend Peterson, A., Papeş, M. & Eaton, M. Transferability and model evaluation in ecological niche modeling: A comparison of GARP and Maxent. Ecography (Cop.). 30, 550–560 (2007).Article 

Google Scholar 
22.Ganeshaiah, K. N. et al. Predicting the potential geographical distribution of the sugarcane woolly aphid Using GARP and DIVA-GIS. Curr. Sci. 85, 1526–1528 (2003).
Google Scholar 
23.Underwood, E. C., Klinger, R. & Moore, P. E. Predicting patterns of non-native plant invasions in Yosemite National Park, California, USA. Divers. Distrib. 10, 447–459 (2004).Article 

Google Scholar 
24.Elith, J. et al. Novel methods improve prediction of species’ distributions from occurrence data. Ecography (Cop.) 29, 129–151 (2006).Article 

Google Scholar 
25.Marmion, M., Parviainen, M., Luoto, M., Heikkinen, R. K. & Thuiller, W. Evaluation of consensus methods in predictive species distribution modelling. Divers. Distrib. 15, 59–69 (2009).Article 

Google Scholar 
26.Araújo, M. B. & New, M. Ensemble forecasting of species distributions. Trends Ecol. Evol. 22, 42–47 (2007).PubMed 
Article 

Google Scholar 
27.Bhatta, M., Shah, K., Devkota, B., Paudel, R. & Panthi, S. Distribution and habitat preference of Red Panda (Ailurus fulgens fulgens) in Jumla District, Nepal. Open J. Ecol. 04, 989–1001 (2014).Article 

Google Scholar 
28.Bista, D. et al. Distribution and habitat use of red panda in the Chitwan-Annapurna Landscape of Nepal. PLoS One 12, e0178797 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
29.Bista, R. & Aryal, A. Status of the Asiatic black bear Ursus thibetanus in the southeastern region of the Annapurna Conservation Area, Nepal. Zool. Ecol. 23 (2013).30.Garshelis, D. & Steinmetz, R. Ursus thibetanus. (errata version published in 2017) The IUCN Red List of Threatened Species. 2016: e. T22824A114252336. (2016).31.Bista, M., Panthi, S. & Weiskopf, S. R. Habitat overlap between Asiatic black bear Ursus thibetanus and red panda Ailurus fulgens in Himalaya Habitat overlap between Asiatic black bear Ursus thibetanus and red panda Ailurus fulgens in Himalaya. https://doi.org/10.1371/journal.pone.0203697 (2018).32.CITES. Asiatic Black bear. Convention on International Trade in Endangered Species of Wild Fauna and Flora https://www.cites.org/eng/gallery/species/mammal/Asiatic_black_bear.html (2019a).33.CITES. Lesser Panda. Convention on International Trade in Endangered Species of Wild Fauna and Flora https://www.cites.org/eng/gallery/species/mammal/lesser_panda.html (2019b).34.Garshelis, Scheick, B., Doan-Crider, D., Beecham & Obbard, M. Ursus americanus, American Black Bear. The IUCN Red List of Threatened Species 2016: e.T41687A45034604. (2016). https://doi.org/10.2305/IUCN.UK.2016-3.RLTS.T41687A45034604.en.35.Chhetri, M. Distribution and abundance of Himalayan black bear and brown bear conflict in Manaslu conservation area. https://ntnc.org.np/index.php/publication/distribution-and-ambundance-himalayan-black-bear-and-brown-bear-and-human-bear-conflict (2013).36.Ali, A. et al. An assessment of food habits and altitudinal distribution of the Asiatic black bear (Ursus thibetanus) in the Western Himalayas, Pakistan. J. Nat. Hist. 51, 689–701 (2017).Article 

Google Scholar 
37.Glatston, A., Wei, F., Zaw, T. & Sherpa, A. P. IUCN red list of threatened species: Ailurus fulgens. (2015).38.Hu, Y. et al. Genomic evidence for two phylogenetic species and long-term population bottlenecks in red pandas. Sci. Adv. 6, eaax5751 (2020).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
39.Chakraborty, R. et al. Status, abundance, and habitat associations of the red panda (Ailurus fulgens) in Pangchen Valley, Arunachal Pradesh, India. Mammalia 79, 25–32 (2015).
Google Scholar 
40.Dorji, S., Vernes, K. & Rajaratnam, R. Habitat correlates of the red panda in the temperate forests of Bhutan. PLoS One 6, e26483 (2011).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
41.Panthi, S., Aryal, A., Raubenheimer, D., Lord, J. & Adhikari, B. Summer diet and distribution of the Red Panda (Ailurus fulgens fulgens) in Dhorpatan hunting reserve, Nepal. Zool. Stud. 51, 701–709 (2012).
Google Scholar 
42.Pradhan, S., Saha, G. K. & Khan, J. A. Ecology of the red panda Ailurus fulgens in the Singhalila National Park, Darjeeling, India. Biol. Conserv. 98, 11–18 (2001).Article 

Google Scholar 
43.Acharya, K. P., Paudel, P. K., Neupane, P. R. & Köhl, M. Human-wildlife conflicts in Nepal: Patterns of human fatalities and injuries caused by large mammals. PLoS One 11, e0161717 (2016).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
44.Liu, Z. et al. Habitat suitability assessment of blue sheep in Helan Mountain based on MAXENT modeling. Acta Ecol. Sin. 33, 7243–7249 (2013).Article 

Google Scholar 
45.Bhusal, N. P. Buffer zone management system in protected areas of Nepal. Third Pole J. Geogr. Educ. 34–44 (2012).46.Carpenter, C. & Zomer, R. Forest ecology of the Makalu-Barun National Park and conservation area, Nepal. Mt. Res. Dev. 16, 135–148 (1996).Article 

Google Scholar 
47.Bhuju, U. R., Shakya, P. R., Basnet, T. B. & Shrestha, S. Nepal biodiversity resource book: Protected areas, Ramsar sites, and World Heritage sites. (International Centre for Integrated Mountain Development (ICIMOD), 2007).48.Wikipedia. Makalu Barun National Park. https://en.wikipedia.org/w/index.php?title=Makalu_Barun_National_Park&oldid=1022613383 (2020).49.Bista, M., Panthi, S. & Weiskopf, S. R. Habitat overlap between Asiatic black bear Ursus thibetanus and red panda Ailurus fulgens in Himalaya. PLoS ONE 13, e0203697 (2018).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
50.Chen, X. & Lei, Y. Effects of sample size on accuracy and stability of species distribution models. A Comparison of GARP and Maxent BT – Recent Advances in Computer Science and Information Engineering, Volume 2. in (eds. Qian, Z. et al.) 601–609 (Springer, Berlin Heidelberg, 2012). https://doi.org/10.1007/978-3-642-25789-6_80.Chapter 

Google Scholar 
51.Zomer, R., Ustin, S. & Ives, J. Using satellite remote sensing for DEM extraction in complex mountainous terrain: Landscape analysis of the Makalu Barun National Park of eastern Nepal. Int. J. Remote Sens. 23, 125–143 (2002).ADS 
Article 

Google Scholar 
52.Shao, Y. & Lunetta, R. S. Comparison of support vector machine, neural network, and CART algorithms for the land-cover classification using limited training data points. ISPRS J. Photogramm. Remote Sens. 70, 78–87 (2012).ADS 
Article 

Google Scholar 
53.Fick, S. E. & Hijmans, R. J. WorldClim 2: New 1-km spatial resolution climate surfaces for global land areas. Int. J. Climatol. 37, 4302–4315 (2017).Article 

Google Scholar 
54.Merow, C., Smith, M. J. & Silander, J. A. Jr. A practical guide to MaxEnt for modeling species’ distributions: What it does, and why inputs and settings matter. Ecography (Cop.) 36, 1058–1069 (2013).Article 

Google Scholar 
55.Steven, J. P., Miroslav, D. & Robert, E. S. Maxent software for modeling species niches and distributions (Version 3.4.1). http://biodiversityinformatics.amnh.org/open_source/Maxent/.56.Phillips, S. J. Transferability, sample selection bias and background data in presence-only modelling: A response to Peterson et al. (2007). Ecography (Cop.) 31, 272–278 (2008).Article 

Google Scholar 
57.Boral, D. & Moktan, S. Predictive distribution modeling of Swertia bimaculata in Darjeeling-Sikkim Eastern Himalaya using MaxEnt: Current and future scenarios. Ecol. Process. 10, 1–16 (2021).Article 

Google Scholar 
58.Pasquale, G. D. et al. Coastal pine-oak glacial refugia in the Mediterranean basin: A biogeographic approach based on charcoal analysis and spatial modelling. Forests 11, 673 (2020).Article 

Google Scholar 
59.Barbet-Massin, M., Jiguet, F., Albert, C. & Thuiller, W. Selecting pseudo-absences for species distribution models: How, where and how many?. Methods Ecol. Evol. 3, 327–338 (2012).Article 

Google Scholar 
60.Adjemian, J. C. Z., Girvetz, E. H., Beckett, L. & Foley, J. E. Analysis of genetic algorithm for rule-set production (GARP) modeling approach for predicting distributions of fleas implicated as vectors of Plague, Yersinia pestis, California. J. Med. Entomol. 43, 93–103 (2006).PubMed 

Google Scholar 
61.Barro, A. S. et al. Redefining the Australian Anthrax Belt: Modeling the Ecological Niche and Predicting the Geographic Distribution of Bacillus anthracis. PLoS Negl. Trop. Dis. 10, e0004689 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
62.Anderson, R. P., Lew, D. & Peterson, A. T. Evaluating predictive models of species’ distributions: Criteria for selecting optimal models. Ecol. Model. 162, 211–232 (2003).Article 

Google Scholar 
63.Babar, S., Giriraj, A., Reddy, C. S., Jentsch, A. & Sudhakar, S. Species distribution models: Ecological explanation and prediction of an endemic and endangered plant species (Pterocarpus santalinus L.f.). Curr. Sci. 102, 1157–1165 (2012).
Google Scholar 
64.Stohlgren, T. J. et al. Ensemble habitat mapping of invasive plant species. Risk Anal. 30, 224–235 (2010).PubMed 
Article 

Google Scholar 
65.Smeraldo, S. et al. Generalists yet different: Distributional responses to climate change may vary in opportunistic bat species sharing similar ecological traits. Mamm. Rev. (2021).66.Pearce, J. & Ferrier, S. Evaluating the predictive performance of habitat models developed using logistic regression. Ecol. Model. 133, 225–245 (2000).Article 

Google Scholar 
67.Chikerema, S., Gwitira, I., Murwira, A., Pfukenyi, D. & Matope, G. Comparison of GARP and Maxent in modelling the geographic distribution of Bacillus anthracis in Zimbabwe. Zimbabwe Vet. J. 35, 1–6 (2017).
Google Scholar 
68.Ray, D., Behera, M. D. & Jacob, J. Evaluating ecological niche models: A comparison between Maxent and GARP for predicting distribution of Hevea brasiliensis in India. Proc. Natl. Acad. Sci. India Sect. B Biol. Sci. 88, 1337–1343 (2018).
Google Scholar 
69.Phillips, S. J. A brief tutorial on Maxent. AT&T Res. 190, 231–259 (2005).
Google Scholar 
70.Jnawali, S. R. et al. The Status of Nepal’s Mammals: The National Red List Series-IUCN (2011).71.Panthi, S., Wang, T., Sun, Y. & Thapa, A. An assessment of human impacts on endangered red pandas (Ailurus fulgens) living in the Himalaya. Ecol. Evol. 9, 13413–13425 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
72.Pearson, R. G. et al. Model-based uncertainty in species range prediction. J. Biogeogr. 33, 1704–1711 (2006).Article 

Google Scholar 
73.Randin, C. F. et al. Are niche-based species distribution models transferable in space?. J. Biogeogr. 33, 1689–1703 (2006).Article 

Google Scholar 
74.Panthi, S., Aryal, A. & Coogan, S. C. P. Diet and macronutrient niche of Asiatic black bear (Ursus thibetanus) in two regions of Nepal during summer and autumn. Ecol. Evol. 9, 3717–3727 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
75.Thapa, A. et al. The endangered red panda in Himalayas: Potential distribution and ecological habitat associates. Glob. Ecol. Conserv. 21, e00890 (2020).Article 

Google Scholar 
76.Shailendra. Human-Bear Conflicts Biological Research Himalayan Black Bear Discovered in Babai Valley of Bardia National. 26, 1999–2001 (2017).77.Acharya, K. P. et al. Pervasive human disturbance on habitats of endangered red panda Ailurus fulgens in the central Himalaya. Glob. Ecol. Conserv. 15, e00420 (2018).Article 

Google Scholar 
78.Letro, L., Wangchuk, S. & Dhendup, T. Distribution of Asiatic black bear and its interaction with humans in Jigme Singye Wangchuck National Park, Bhutan. Nat. Conserv. Res. 5, 44–52 (2020).Article 

Google Scholar 
79.Karki, S. T. Do protected areas and conservation incentives contribute to sustainable livelihoods? A case study of Bardia National Park, Nepal. 988–999.80.Guisan, A. & Zimmermann, N. E. Predictive habitat distribution models in ecology. Ecol. Model. 135, 147–186 (2000).Article 

Google Scholar  More

SoilContents of risk metals in soilsLands of localities from which soil and plant samples were taken belong to agricultural lands.Soil reaction is one of the factors that most affects the behaviour of heavy metals in soil. Low pH values pose a risk of reduced nutrient intake and increase the availability of heavy metals for plants29,30.The presence of risk elements in the soil was evaluated based on their contents in bioavailable form (mobile forms), determined in soil extracts NH4NO3, and the total contents of risk elements were determined in soil extract by aqua regia (Table 1).Table 1 The contents of risk elements (Cu, Ni, Pb, Cd, Hg, Mn) in soil (mg/kg).Full size tableAccessible heavy metals for plants are those which are present in the soil solution as soluble components or those which are easily dissolved by root exudates31. The highest Cu contents determined in soil extract by NH4NO3, were in the cadastre of Gemerská Poloma (max. 0.390 mg/kg) (Table 1). However, even the highest determined concentration of Cu in its bioavailable form did not exceeded the determined critical value for this element18. Nickel is a beneficial element for plants. Elevated Ni concentrations in soils have a potential negative effect on plants32. Content of bioavailable forms of nickel is lower than the determined critical value in all analysed samples. Cadmium and lead present a risk to agricultural activity in this area. Cadmium in soil is highly bioavailable and has higher mobility in plants compared to other heavy metals. It is easily transported by roots to shoots. In contrast, lead is one of the least mobile heavy metals. It is naturally concentrated in the upper layers of the soil33. The contents of the available forms of cadmium and lead exceed the critical values for these elements. In case of lead, the determined contents are from 0.257 Henckovce to 0.676 Gemerská Poloma. Takáč et al.34 determined in 20 soil samples from the Central Spiš region 7.2–257.6 mg Cu/kg soil and 1.0–84.8 mg Pb/kg in their potentially mobilizable form and 0.4–1.4 mg Cu/kg soil and 4.3–7.1 mg Pb/kg in their mobile form. In comparison with our results, Vilček et al.35 determined a lower content of Cd (0.04), Pb (0.17), Ni (0.15) and higher Cu content (0.48) mg/kg in forms accessible to plants in 16 soil samples from locality Nižná Slaná in the years 2006–2008. However, high concentrations of metals in soil do not necessarily mean the availability of metals for plants36. As a result, extractable Mn is often a better indicator of Mn availability. Mn2+ is generally considered to be bioavailable22. The highest concentration of Mn was measured in soil samples from the cadastre of Nižná Slaná. On the contrary, the lowest concentrations were detected in samples from Gemerská Poloma cadastre, which is the furthest cadastre from the source. No critical limit is set up for manganese according to Slovak legislation, it is not possible to classify these soils as contaminated/uncontaminated. For comparison, the EDTA-extractable content of Mn ranged from 22.7 to 127 mg/kg dry soil (China)29; the mobile concentrations between 0.32 and 202.0 mg/kg and the available concentrations from 5.4 to 126.3 mg/kg (Egypt)37.Based on results of statistical analysis, significant higher content of Cu, Pb and Cd can be stated in samples from Gemerská Poloma cadastre. These soils are classified as gley fluvisols, soils from the other two localities are cambisols (from medium heavy to light) and acid cambisols (Henckovce), cambisols from medium heavy to light and typically acid cambisols (Nižná Slaná). The soil profile of fluvisols is repeatedly disrupted by floods, which often enriches them with a new layer of sludge sediments2.Another method for determination of metal content in soil is mineralisation using aqua regia, which dissolves most of the soil constituents except those strongly bound in silicate minerals. This content is sometimes referred to as pseudototal (determined in aqua regia). In this way, all elements that are likely to become bioavailable in the long term are determined38.Pseudototal contents of risk metals (Table 1) determined in soil extract using aqua regia were higher than their limit value in case of Cu (Gemerská Poloma cadastre), Cd (all cadastres) and Hg (cadastre of Henckovce and Gemerská Poloma).Due to the fact that the hygienic condition of agricultural soils is assessed according to the exceeding of the limit values of at least one risk substance, the monitored plots can be classified as contaminated (Cu  > 60.0, Cd  > 0.7, Hg  > 0.5 mg/kg soil).Manganese is not classified as risk element in Slovak legislation.Tóth et al.39 classified European soils into four categories: (1) no detectable content of HM, (2) the concentration of the investigated element is above the threshold value (Hg 0.5, Cd 1, Cu 100, Pb 60 and Ni 50 mg/kg), but below the lower guideline value (Hg 2, Cd 10, Cu 150, Pb 200 and Ni 100 mg/kg), (3) concentration of one or more element exceeds the lower guideline value but is below the higher guideline value (Hg 5, Cd 20, Cu 200, Pb 750 and Ni 150 mg/kg), (4) samples having concentrations above the higher guideline value.In comparison with the threshold and guideline values, soils in cadastres of Gemerská Poloma (Cu), Henckovce, Nižná Slaná, Gemerská Poloma (Cd, Hg) represent the ecological risk. Threshold and guideline values for Mn were not defined.The Spiš region and the northern part of the Gemer region belong to the most polluted areas in Slovakia in terms of soil contamination due to mining and metallurgical activities that have been carried out here in the past. Soils near the sludge in Nižná Slaná contain 3.17–53.3 (14.2–301, 0.71–20.6, 3.33–177, 12.9–223 and 675–11,510, respectively) mg Cd (Cu, Hg, Ni, Pb and Mn, respectively)/kg of soil14. In loaded area of Dongchuan, (China), contained Cd (Cu, Hg, Ni and Pb, resp.) 0.20–3.57 (45.38–2026, 0.02–0.23, 24.06–95.9 and 6.83–146.6, resp.) mg/kg40. In contrast, in the agricultural area of Punjab of the India, the soil contamination was caused by an excessive use of agrochemicals and polluted irrigation sources. Increased Cu (Pb and Cd) contents were determined in the soil samples: 9.0–48.5 (5.5–9.67 and 0.516–1.58, resp.) mg/kg41.However, in most cases, a large portion of the total element content is not available for immediate uptake by plants. Available forms represent a small proportion of this total content which is potentially available to plants. Availability is affected by many factors, including pH, redox state, macronutrient levels, available water content and temperature29,33,36,38.Indicators of soil contaminationContamination factors and degree of contaminationThe contamination character may be described in a uniform, adequate and standardised way by means of the contamination factor and the degree of contamination. Hakanson24 reported four Contamination degrees of individual metal (({mathrm{C}}_{mathrm{f}}^{mathrm{i}})) – low (({mathrm{C}}_{mathrm{f}}^{mathrm{i}}) < 1), moderate (1 ≤ ({mathrm{C}}_{mathrm{f}}^{mathrm{i}})   More