Ecology

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

1.Madsen, E. L. Identifying microorganisms responsible for ecologically significant biogeochemical processes. Nat. Rev. Micro. 3, 439 (2005).CAS 
Article 

Google Scholar 
2.Bell, T., Newman, J. A., Silverman, B. W., Turner, S. L. & Lilley, A. K. The contribution of species richness and composition to bacterial services. Nature 436, 1157–1160 (2005).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
3.Delgado-Baquerizo, M. et al. Microbial diversity drives multifunctionality in terrestrial ecosystems. Nat. Commun. 7, 10541 (2016).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
4.Galand, P. E., Pereira, O., Hochart, C., Auguet, J. C. & Debroas, D. A strong link between marine microbial community composition and function challenges the idea of functional redundancy. ISME J. 12, 2470 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
5.Galand, P. E., Salter, I. & Kalenitchenko, D. Ecosystem productivity is associated with bacterial phylogenetic distance in surface marine waters. Mol. Ecol. 24, 5785–5795 (2015).PubMed 
Article 
PubMed Central 

Google Scholar 
6.Chase, J. M. Community assembly: When should history matter?. Oecologia 136, 489–498 (2003).ADS 
PubMed 
Article 
PubMed Central 

Google Scholar 
7.Lozupone, C. A. & Knight, R. Global patterns in bacterial diversity. Proc. Natl. Acad. Sci. USA 104, 11436–11440 (2007).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
8.Thompson, L. R. et al. A communal catalogue reveals Earth’s multiscale microbial diversity. Nature 551, 457 (2017).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
9.Hanson, C. A., Fuhrman, J. A., Horner-Devine, M. C. & Martiny, J. B. Beyond biogeographic patterns: Processes shaping the microbial landscape. Nat. Rev. Micro. 10, 497 (2012).CAS 
Article 

Google Scholar 
10.Sunagawa, S. et al. Structure and function of the global ocean microbiome. Science 348, 1261359 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
11.Fukami, T. Historical contingency in community assembly: Integrating niches, species pools, and priority effects. Annu. Rev. Ecol. Evol. Syst. 46, 1–23 (2015).Article 

Google Scholar 
12.Martiny, J. B. H. et al. Microbial biogeography: Putting microorganisms on the map. Nat. Rev. Micro. 4, 102 (2006).CAS 
Article 

Google Scholar 
13.Hawkes, C. V. & Keitt, T. H. Resilience vs. historical contingency in microbial responses to environmental change. Ecol. Lett. 18, 612–625 (2015).PubMed 
Article 
PubMed Central 

Google Scholar 
14.Bouskill, N. J. et al. Pre-exposure to drought increases the resistance of tropical forest soil bacterial communities to extended drought. ISME J. 7, 384 (2013).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
15.Ge, Y. et al. Differences in soil bacterial diversity: Driven by contemporary disturbances or historical contingencies?. ISME J. 2, 254 (2008).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
16.Rousk, J., Smith, A. R. & Jones, D. L. Investigating the long-term legacy of drought and warming on the soil microbial community across five European shrubland ecosystems. Glob. Change Biol. 19, 3872–3884 (2013).ADS 
Article 

Google Scholar 
17.Langenheder, S., Lindström, E. S. & Tranvik, L. J. Structure and function of bacterial communities emerging from different sources under identical conditions. Appl. Environ. Microbiol. 72, 212–220 (2006).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
18.Langenheder, S., Lindström, E. S. & Tranvik, L. J. Weak coupling between community composition and functioning of aquatic bacteria. Limnol. Oceanogr. 50, 957–967 (2005).ADS 
Article 

Google Scholar 
19.Vass, M. & Langenheder, S. The legacy of the past: Effects of historical processes on microbial metacommunities. Aquat. Microb. Ecol. 79, 13–19 (2017).Article 

Google Scholar 
20.Svoboda, P., Lindström, E. S., Osman, O. A. & Langenheder, S. Dispersal timing determines the importance of priority effects in bacterial communities. ISME J. 12, 644 (2018).PubMed 
Article 
PubMed Central 

Google Scholar 
21.Rummens, K., De Meester, L. & Souffreau, C. Inoculation history affects community composition in experimental freshwater bacterioplankton communities. Environ. Microbiol. 20, 1120–1133 (2018).PubMed 
Article 
PubMed Central 

Google Scholar 
22.Andersson, M. G., Berga, M., Lindström, E. S. & Langenheder, S. The spatial structure of bacterial communities is influenced by historical environmental conditions. Ecology 95, 1134–1140 (2014).PubMed 
Article 
PubMed Central 

Google Scholar 
23.Vyverman, W. et al. Historical processes constrain patterns in global diatom diversity. Ecology 88, 1924–1931 (2007).PubMed 
Article 
PubMed Central 

Google Scholar 
24.Sefbom, J., Sassenhagen, I., Rengefors, K. & Godhe, A. Priority effects in a planktonic bloom-forming marine diatom. Biol. Lett. 11, 20150184 (2015).PubMed 
PubMed Central 
Article 

Google Scholar 
25.Kalenitchenko, D. et al. Ecological succession leads to chemosynthesis in mats colonizing wood in sea water. ISME J. 10, 2246–2258 (2016).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
26.Kalenitchenko, D., Le Bris, N., Peru, E. & Galand, P. E. Ultrarare marine microbes contribute to key sulphur-related ecosystem functions. Mol. Ecol. 27, 1494–1504 (2018).PubMed 
Article 
PubMed Central 

Google Scholar 
27.Ghiglione, J. F. et al. Role of environmental factors for the vertical distribution (0–1000 m) of marine bacterial communities in the NW Mediterranean Sea. Biogeosciences 5, 1751–1764 (2008).ADS 
CAS 
Article 

Google Scholar 
28.Ghiglione, J.-F. et al. Pole-to-pole biogeography of surface and deep marine bacterial communities. Proc. Natl. Acad. Sci. USA 109, 17633–17638 (2012).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
29.Salazar, G. et al. Gene expression changes and community turnover differentially shape the global ocean metatranscriptome. Cell 179, 1068-1083.e1021 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
30.Kalenitchenko, D. et al. The early conversion of deep-sea wood falls into chemosynthetic hotspots revealed by in situ monitoring. Sci. Rep. 8, 907. https://doi.org/10.1038/s41598-017-17463-2 (2018).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
31.Kalenitchenko, D. et al. Temporal and spatial constraints on community assembly during microbial colonization of wood in seawater. ISME J. 9, 2657–2670 (2015).PubMed 
PubMed Central 
Article 

Google Scholar 
32.Kalenitchenko, D. et al. Bacteria alone establish the chemical basis of the wood-fall chemosynthetic ecosystem in the deep-sea. ISME J. 12, 367–379 (2018).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
33.Galand, P., Salter, I. & Kalenitchenko, D. Microbial productivity is associated with phylogenetic distance in surface marine waters. Mol. Ecol. 24, 5785–5795 (2015).PubMed 
Article 
PubMed Central 

Google Scholar 
34.Edgar, R. C. Search and clustering orders of magnitude faster than BLAST. Bioinformatics 26, 2460–2461. https://doi.org/10.1093/bioinformatics/btq461 (2010).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
35.Pruesse, E. et al. SILVA: a comprehensive online resource for quality checked and aligned ribosomal RNA sequence data compatible with ARB. Nucleic Acids Res. 35, 7188–7196. https://doi.org/10.1093/nar/gkm864 (2007).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
36.Cox, M. P., Peterson, D. A. & Biggs, P. J. SolexaQA: At-a-glance quality assessment of Illumina second-generation sequencing data. BMC Bioinform. 11, 485. https://doi.org/10.1186/1471-2105-11-485 (2010).Article 

Google Scholar 
37.Rho, M., Tang, H. & Ye, Y. FragGeneScan: Predicting genes in short and error-prone reads. Nucleic Acids Res. https://doi.org/10.1093/nar/gkq747 (2010).Article 
PubMed 
PubMed Central 

Google Scholar 
38.Wilke, A. et al. The M5nr: A novel non-redundant database containing protein sequences and annotations from multiple sources and associated tools. BMC Bioinform. 13, 141. https://doi.org/10.1186/1471-2105-13-141 (2012).CAS 
Article 

Google Scholar 
39.Kanehisa, M., Sato, Y., Kawashima, M., Furumichi, M. & Tanabe, M. KEGG as a reference resource for gene and protein annotation. Nucleic Acids Res 44, D457–D462 (2016).CAS 
Article 

Google Scholar 
40.Anders, S. & Huber, W. Differential expression analysis for sequence count data. Nat Précéd 1–1 https://doi.org/10.1038/npre.2010.4282.1 (2010).Article 

Google Scholar 
41.Meyer, F. et al. The metagenomics RAST server—A public resource for the automatic phylogenetic and functional analysis of metagenomes. BMC Bioinform. 9, 386. https://doi.org/10.1186/1471-2105-9-386 (2008).CAS 
Article 

Google Scholar 
42.Dixon, P. VEGAN, a package of R functions for community ecology. J Veg Sci 14, 927–930 (2003).Article 

Google Scholar 
43.Blanchette, R. A., Nilsson, T., Daniel, G. & Abad, A. Biological Degradation of Wood. in vol. 225, 141–174 (American Chemical Society, 1989).
Google Scholar 
44.Fagervold, S. K. et al. Microbial communities associated with the degradation of oak wood in the Blanes submarine canyon and its adjacent open slope (NW Mediterranean). Prog. Oceanogr. 118, 137–143. https://doi.org/10.1016/j.pocean.2013.07.012 (2013).ADS 
Article 

Google Scholar 
45.Sommer, U. Convergent succession of phytoplankton in microcosms with different inoculum species composition. Oecologia 87, 171–179 (1991).ADS 
PubMed 
Article 
PubMed Central 

Google Scholar 
46.Weiher, E. & Keddy, P. A. The assembly of experimental wetland plant communities. Oikos 73, 323–335 (1995).Article 

Google Scholar 
47.Wilson, J. B. et al. A test of community reassembly using the exotic communities of New Zealand roadsides in comparison to British roadsides. J. Ecol. 88, 757–764 (2000).Article 

Google Scholar 
48.Kodric-Brown, A. & Brown, J. H. Highly structured fish communities in Australian desert springs. Ecology 74, 1847–1855 (1993).Article 

Google Scholar 
49.Grover, J. P. & Lawton, J. H. Experimental studies on community convergence and alternative stable states: Comments on a paper by Drake et al. J. Anim. Ecol. 63, 484–487 (1994).Article 

Google Scholar 
50.Lawler, S. P. Direct and indirect effects in microcosm communities of protists. Oecologia 93, 184–190 (1993).ADS 
PubMed 
Article 
PubMed Central 

Google Scholar 
51.Chase, J. M. Experimental evidence for alternative stable equilibria in a benthic pond food web. Ecol. Lett. 6, 733–741 (2003).Article 

Google Scholar 
52.Petraitis, P. S. & Latham, R. E. The importance of scale in testing the origins of alternative community states. Ecology 80, 429–442 (1999).Article 

Google Scholar 
53.Hiscox, J. et al. Priority effects during fungal community establishment in beech wood. ISME J. 9, 2246 (2015).PubMed 
PubMed Central 
Article 

Google Scholar 
54.Fukami, T. et al. Assembly history dictates ecosystem functioning: evidence from wood decomposer communities. Ecol. Lett. 13, 675–684 (2010).PubMed 
Article 
PubMed Central 

Google Scholar 
55.Dhami, M. K., Hartwig, T. & Fukami, T. Genetic basis of priority effects: Insights from nectar yeast. Proc. R. Soc. Lond. B. 283, 20161455 (2016).
Google Scholar 
56.Fukami, T. & Morin, P. J. Productivity–biodiversity relationships depend on the history of community assembly. Nature 424, 423 (2003).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
57.Khelaifia, S. et al. Desulfovibrio piezophilus sp. nov., a piezophilic, sulfate-reducing bacterium isolated from wood falls in the Mediterranean Sea. Int. J. Syst. Evol. Micr. 61, 2706–2711 (2011).CAS 
Article 

Google Scholar 
58.Sievert, S. M., Wieringa, E. B., Wirsen, C. O. & Taylor, C. D. Growth and mechanism of filamentous-sulfur formation by Candidatus Arcobacter sulfidicus in opposing oxygen-sulfide gradients. Environ. Microbiol. 9, 271–276 (2007).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar  More

A simple methodological framework was established for testing the BCR model using empirical datasets, consisting of the GPS data of 5 animals. For each of these 5 animals, the three parameters were accordingly tuned using a straightforward estimation procedure. This procedure uses the empirical datasets to infer the parameters’values (Fig. 1). We also used the datasets to assess the model’s reliability—or performance-. We also detail other analyzes that were carried out to ensure the robustness and consistency of the approach, including the deterministic nature of the 5 statistics and a sensitivity analysis. This analysis consists in evaluating the performance of the BCR using a sweep method that produce arbitrary values of the parameters instead of using data-driven estimations. All BCR simulations and the five statistics were performed using MATLAB Version 7.13.0.564 (R2011b).Figure 1Framework used for testing the BCR model performance, for one animal. Black lines detail the two operations processed from the GPS dataset. The 3 parameters are estimated from the GPS data and—using these parameters—1000 simulations of the BCR model are computed. No particular operations are associated with the dotted black lines, but they show how the BCR and the GPS dataset are evaluated and compared using the statistics.Full size imageDataThe locations of 5 GPS-collared red deer (Cervus elaphus) were gathered at La Petite Pierre National Hunting and Wildlife Reserve (NHWR), in north-eastern of France (48.8321 (Lat.) / 7.3514 (Lon.)). The reserve is an unfenced 2670 ha forest area characteristics by deciduous trees (mostly Fagus sylvatica) in the western part and by coniferous species (mostly Pinus sylvestris and Abies alba) in the eastern part in nature reserve surrounded by crops and pastures. It is located at a low elevation area of the Vosges mountain range, which rises up to 400 m a. s. l. The climate is continental with cool summers and mild winters (mean January and July temperatures of 1.4 and 19.6 (^{circ })C, respectively, data from Phalsbourg weather station, Meteo France, from 2004 to 2017). Three ungulate species are present and mainly managed through hunting in the NHWR: wild boar, red deer and roe deer. The present study focuses on female red deer for test model. A detailed overview of the landscape and surroundings is given in40. The GPS data had regular observation frequencies with high frequency sampling (Table 1). In the following text, we note (X_i = [X_i^{(1)}, X_i^{(2)}]) the locations of the individual with (X_i in {mathbb {R}}^{2}), (i=1,2,ldots ,n) and where (X_i^{(1)}), (X_i^{(2)}) represent the longitude and latitude respectively. We use (t_{i}) ((t_1=0)) as the time elapsed between two successive locations (X_{i-1}), (X_{i}) and$$begin{aligned} {overline{T}}= dfrac{1}{n} sum _{i=1}^{n} t_i end{aligned}$$
(1)
as the average sampling time. The trajectory of the animal, or ‘path’, was interpolated using linear interpolation between each pair of recorded observations (Fig. 2 and detailed in Supplementary Methods (Eq. 21) and associated Graphic 2). It approximates the animal travels in straight lines at constant velocity between each pair of locations41. The attractor (X_F) of one individual was estimated as the isobarycenter of all recorded locations:$$begin{aligned} X_F = left[ frac{1}{n}sum _{i=1}^{n} X_i^{(1)}, frac{1}{n}sum _{i=1}^{n} X_i^{(2)} right] end{aligned}$$
(2)
Table 1 Data summary. For each animal, the total number of observations n is given along with the period of collection (date and time), the sampling rate ({overline{T}}) (i.e. the average time between 2 observations) (in min.) and corresponding standard deviation, total distance (in kilometers), total recording time (in days) and average speed s (in (10^{-2}) m.s-1).Full size tableFigure 2Individual paths of the five red deer. Individual paths of the five red deer. The individual paths are plotted for the five red deer (left panel, a) along with the distribution of the relative turning angles (degrees) in polar plots (right panel, b). An angular value of 0 consists in a straight motion from the previous location, while a relative turning angle of 180 (^{circ }) c corresponds to a turn back.Full size image
BCR modelThe model aims at estimating the location at the next time step, given the actual location X at step i:$$begin{aligned} X_{i+1} = f left( X_i right) end{aligned}$$
(3)
such that the function (f(cdot )) is assumed to be representative of the behavior of the animal on sufficiently large time scales. We considered one individual of a given species with no interaction and simulated its movement in continuous space and discrete time in 2 dimensions. The BCR includes 3 parameters coupled with isotropic diffusion:

Diffusion: A random direction with uniform spatial distribution in a 2D plane,

Bias ((p_F)): An increased probability to go to a fixed point named attractor42. This attractor was estimated as the isobarycenter of all recorded locations, defined as (X_F) (Eq. 2). This yields a bias or advection parameter in the direction of (X_F). We use the term ’attraction’ for the bias component of the BCR and the term ’den’ for the attractor. In the data set we study, the den is equivalent to the deer’s bunk.

Correlated component ((p_I)): This parameter increases the probability to move forward, i.e. to perform one step in the direction of the previous step. This is equivalent to a short term bias in movement, when the animal has inertia. We refers to ’inertia’ for the correlated component,

Immobility ((p_s)): We included this as a specific parameter and the movement is stopped for one step. This takes into account the absence of movement between a pair of locations. It can be accredited to technological limitations with the satellite telemetry due to a weak GPS signal strength, possibly due to natural elements: such as when the animal was standing underneath a rock or due to dense clouds, dust particles, mountains or flying objects, such as airplanes). However, this can also be part of the behavior of the animals, during specific times: sleep cycles or foraging for instance. We use (d_{min }) to denote this distance cutoff and set (d_{min }=10)m which corresponds to the magnitude of the error typically found in GPS locations43. We also use (d_{min }) to encapsulate GPS error and peculiar ecological behavior, not associated with (p_I) or (p_F), that are beyond the scope of this study.

The effect of each parameter is detailed in Fig. 3. The typical model contains all three parameters: (p_I), (p_s) and (p_F) for describing animal motion while offering a trade-off between the number of parameters and the description of animal motion.Figure 3Simulated animal motions over arbitrary parameter values. Fifty motions of length (n_s=100) steps are simulated and originate from a common centroid (downward-pointing triangle) with increased levels of correlation ((p_I)), immobility ((p_s)) and bias ((p_F)). Both the location of the attractor (X_F) (black dot) and the log-normal parameters controlling the step size distribution are fixed ((mu =3), (sigma ^2 = 1)).Full size imageEstimation of the parametersThe estimation of the three parameters for each animal is based on the empirical datasets. We distinguished between the states, where one state is described by the pair (left{ X_{i-1}, X_iright} ) and the situations, where one situation is described by the past ((X_{i-1})), current ((X_i)) and future ((X_{i+1})) locations. Knowing both the state of the animal at a given time step i and its situation—the realization of movement at the next time step (i+1)—allowed for collecting the occurrences of inertia, immobilism and attraction. This could be done provided we account for the variability of the movement: the animal may not be heading exactly toward the den, or performing inertia with an exact angular value of (pi ). Thus we discretized the space around the animal in 8 quadrants at each time step i. For example, if the animal was heading straight forward with a margin of (pm pi /8) then it was considered in the situation of inertia. In other words, the state could fall in a situation of inertia with a margin of (pm pi /8). Such a discretization can be represented as a matrix, depending on the state of the animal, its location and the location of the den at each time step (see Supplementary methods, Eq. 19). In order to gather enough data samples per situation, we arbitrary used angular thresholds of (pi /8) as a convenient trade-off between data scarcity and precision loss. Using smaller threshold values (say (pi /10)) may result in too few samples per situations. Using larger threshold values such as (pi /4) may result in a loss of precision while capturing additional movement samples that may not correspond to the situation.We first needed to define in which state is the animal at each step i. A state is the 2-tuple containing the previous and actual observation ({X_{i-1}, X_i}). We wanted to distinguish between non-conflicting and conflicting states, where a non-conflicting state is when the animal is in one state only, while a conflicting state is when the animal is in two states at once. We defined two conflicting states:$$begin{aligned} {mathscr {H}}_{IF}:={i : widehat{ left| X_{i-1} X_i X_F right| } le pi /8} end{aligned}$$
(4)
when the animal was already heading toward the den (X_F), and:$$begin{aligned} {mathscr {H}}_{Is}:={i : dleft( X_{i-1}, X_iright) le d_{min }} end{aligned}$$
(5)
when the distance between two consecutive observations was too small ((le d_{min }) m.), describing an individual that was already immobile. Such that the subset of non-conflicting states is:$$begin{aligned} {mathscr {H}}:={1,cdots ,n} – {mathscr {H}}_{IF} – {mathscr {H}}_{Is} end{aligned}$$
(6)
We then needed to assess in which situation the animal was for each corresponding state. A situation is the 3-tuple (left{ X_{i-1}, X_i, X_{i+1}right} ). We defined three subsets of situations corresponding to a straight forward motion (I), no motion (s) and a motion toward the den (F):$$begin{aligned} I:= {i: pi – pi /8 < widehat{(X_{i-1} X_i X_{i+1})} le pi + pi /8} end{aligned}$$ (7) $$begin{aligned} s:= {i: dleft( X_i, X_{i+1} right) le d_{min }} end{aligned}$$ (8) $$begin{aligned} F:= {i: mid widehat{X_{i-1} X_t X_{i+1}} - widehat{X_{i-1} X_i X_F} mid le pi /8} end{aligned}$$ (9) With (d(cdot ,cdot )) the Euclidean distance between two locations. For the situations in s, we considered that the animal is not performing a motion if the Euclidean distance between two successive locations was (le d_{min })m.We counted the number of states falling in each situation, for states in ({mathscr {H}}) (Eq. 6). We defined (x_{1}), (x_{2}), (x_{3}) as the empirical proportion of cases corresponding to each situation:$$begin{aligned} {left{ begin{array}{ll} x_1 = dfrac{# I cap {mathscr {H}}}{# {mathscr {H}}} ; qquad x_1:=dfrac{1+p_I}{chi } \[16pt] x_2 =dfrac{# s cap {mathscr {H}}}{# {mathscr {H}}} ; qquad x_2:=dfrac{p_s}{chi }\[16pt] x_3 = dfrac{# F cap {mathscr {H}}}{# {mathscr {H}}} ; qquad x_3:=dfrac{1+p_F}{chi } end{array}right. } end{aligned}$$ (10) with (chi = 8+p_{I}+p_{s}+p_{F}). The values of (x_1), (x_2) and (x_3) were then gathered for each animal. We did not use immobile locations (i.e. distances separating two successive observations must be ( > d_{min }) m) for the estimations of (x_1) and (x_3). Solving Eq. (10) for (chi ) with respect to (x_1), (x_2), (x_3) yields:$$begin{aligned} chi = dfrac{6}{1-(x_{1} + x_{2} + x_{3})} end{aligned}$$
(11)
Plugging in Eq. 10:$$begin{aligned} {left{ begin{array}{ll} p_{I} = x_1 chi -1\ p_{s} = x_2 chi \ p_{F} = x_3 chi -1 end{array}right. } end{aligned}$$
(12)
Note that we assumed that (p_{IF} = p_I + p_F) in ({mathscr {H}}_{IF}) and (p_{Is} = p_I + p_s) in ({mathscr {H}}_{Is}) as a convenient arrangement and ignoring higher order conflicting cases. Investigating the step-size distribution in the 5 deers, we found a log-normal step size distribution (Supplementary Fig. S1). We then set a log-normal distribution (ln {mathscr {N}}(mu , sigma ^2)) for the step size distribution for the step size in the BCR.The same estimation procedure was used for configurations using a different number of parameters and quantity (chi ) is accordingly calculated depending on the number of parameters used. It is possible to obtain negative values using this inference method. A parameter with a negative values reflects a direction that is not favored by the animal. In such a case, one should rethink the design of the BCR by changing the parameters (see Supplementary methods, section “negative parameters”). In the subsequent sections, we only consider parameters with positives values.BCR dynamicsThe BCR dynamics for each animal are completely determined by the three parameters (p_F), (p_I), (p_s), taking values in ({mathbb {R}}^{+}), and the step-size distribution. If (p_F = p_I = p_s = 0), the BCR resumes to a typical two-dimensional random walk with a log-normal step size distribution (ln {mathscr {N}}(mu , sigma ^2)). The dynamics can be visualized in Fig. 3 for different values of each parameter. When simulating a step in the model, the motion in ({mathscr {H}}) is described by:$$begin{aligned} f left( X_i right) = {left{ begin{array}{ll} left{ X_i^{(1)} + d cos (alpha _1) ; X_i^{(2)} + d sin (alpha _1) right} &{} qquad text {if } x in [0,8[ \ left{ X_i^{(1)} + d cos (alpha _2) ; X_i^{(2)} + d sin (alpha _2) right} &{} qquad text {if } x in [8,8+p_I[ \ X_i &{} qquad text {if } x in [8+p_I, 8+p_I+p_s[ \ left{ X_i^{(1)} + d cos (alpha _3) ; X_i^{(2)} + d sin (alpha _3) right} &{} qquad text {else} end{array}right. } end{aligned}$$
(13)
with x, d, (alpha _1) random variables defined as (x sim {mathscr {U}} in [0, chi ]), (d sim ln {mathscr {N}}(mu , sigma )), (alpha _1 sim {mathscr {U}} in [0,2pi ]). Variables (alpha _2), (alpha _3) are related to the angular values (alpha _2 = {{,{mathrm{atan2}},}}(X_{i}^2 – X_{i-1}^2, X_{i}^1 – X_{i-1}^1)), (alpha _3 = {{,{mathrm{atan2}},}}(X_{F}^{(2)} – X_{i}^{(2)}, X_{F}^{(1)} – X_{i}^{(1)})) with ({{,{mathrm{atan2}},}}(y, x)) the four quadrant inverse tangent function (14):$$begin{aligned} {{,{mathrm{atan2}},}}(y, x) = {left{ begin{array}{ll} arctan left( {frac{y}{x}}right) &{} x > 0,\ arctan left( {frac{y}{x}}right) +pi &{} x< 0{text {, }}y ge 0,\ arctan left( {frac{y}{x}}right) -pi &{} x< 0{text {, }}y< 0,\ +{frac{pi }{2}} &{} x=0{text {, }}y > 0,\ -{frac{pi }{2}} &{} x=0{text {, }}y < 0,\ 0 &{} x=0{text {, }}y=0text {.} end{array}right. } end{aligned}$$ (14) The motion in ({mathscr {H}}_{Is}) is:$$begin{aligned} f left( X_i right) = {left{ begin{array}{ll} left{ X_i^{(1)} + d cos (alpha _1) ; X_i^{(2)} + d sin (alpha _1) right} &{} qquad text {if } x in [0,8[ \ X_i &{} qquad text {if } x in [8, 8+p_I+p_s[ \ left{ X_i^{(1)} + d cos (alpha _3) ; X_i^{(2)} + d sin (alpha _3) right} &{} qquad text {else} end{array}right. } end{aligned}$$ (15) The motion in ({mathscr {H}}_{IF}) is:$$begin{aligned} f left( X_t right) = {left{ begin{array}{ll} left{ X_i^{(1)} + d cos (alpha _1) ; X_i^{(2)} + d sin (alpha _1) right} &{} qquad text {if } x in [0,8[ \ X_t &{} qquad text {if } x in [8, 8+p_s[ \ left{ X_i^{(1)} + d cos (alpha _2) ; X_i^{(2)} + d sin (alpha _2) right} &{} qquad text {else} end{array}right. } end{aligned}$$ (16) Statistics for describing animal movementWe simulated (N=1000) BCR and used 5 statistics to assess the model reliability on spatial features including: (i) the distribution of relative turning angles which provides information about the movement of the animal, (ii) the home range which provides information about the spatial density of observations and (iii) observation counts using still and mobile transects, providing information on absolute observation abundance44. A detailed description of each statistic is provided in Supplementary Methods and Fig. S2. The reliability—or performance—was assessed in each animal and studied statistic using two error terms (e_1) and (e_2). Error (e_1) is the (L^1) norm to compare the differences between the statistic ({tilde{mathscr {S}}}) computed over a simulated path, and the statistic (smash {{mathscr {S}}}) computed over the data-set:$$begin{aligned} e_1 mathrel {mathop :}= sum {text {errors}} = sum _{k=1}^N |smash {{mathscr {S}}}- {tilde{mathscr {S}}}_k| end{aligned}$$ (17) With (k = 1,ldots , N) the number of simulations of the BCR. Error (e_1) is the sum of absolute differences in the given statistic, and is a natural way of measuring the distance between the statistics computed on the data set and the trajectories generated using the BCR. We also focused on the average relative error (e_2) as an indicator of the sensitivity:$$begin{aligned} e_2 mathrel {mathop :}= dfrac{1}{N} sum _{k=1}^N dfrac{{tilde{mathscr {S}}}_k}{smash {{mathscr {S}}}} end{aligned}$$ (18) Distribution of turning anglesFor each individual, the distribution of counter-clockwise relative turning angles (widehat{(X_{i-1} X_i X_{i+1})}) was gathered, provided (d(X_{i-1}, X_{i}) > d_{min }) and (d(X_{i}, X_{i+1}) > d_{min }). This means that we only kept the angles from observations that were separated by an Euclidean distance greater than (d_{min }).Home rangeWe used an adaptive kernel density estimator (matlab package kde2d—kernel density estimation version 1.3.0.0) as an estimator of the utilization distribution45 to represent the home range of the animal. The approach of Z.I. Botev provided an estimate of observation density using a bivariate (Gaussian) kernel with diagonal bandwidth matrix46. The density was estimated over a grid of (210 times 210) nodes and we computed the home range area (in m2) for various values: 100, 99, 95, 90, 80, …, 20, 10% of the estimated density. Similarly to the distribution of turning angles, we compared each value of the data’s home range against the simulated one.DilationDilation is generally used to account for the spatial attributes of an object such as to measure an area around the path or the volume of a brownian motion (see Wiener sausage47 and Gromov–Hausdorff distance). In our approach, we use dilation of both simulated and GPS paths for two reasons: to have a real—and comparable—number that accounts for how a trajectory has explored space and because it is natural tool from a census point of view (the dilated path corresponds to the area where the animal can be detected). Each simulated or real path was plotted in binary format in a window and dilated with a disk shape. The window size was set to a huge value in order to encapsulate the dilated path while preventing boundary effects, i.e. the convex envelope of the dilated area did not collide with any window border. We then estimated the surface covered by the dilated path for 100 different sizes of the disk, from disk size 1 to disk size 100. We compared each value of the data’s estimated surface against the simulated one.Immobile transectsWe used still transects that counted the number of times the animal was seen in their line of sight. We arbitrary set the line-of-sight value at 200 m. The number of sightings of each transect was gathered and ordered in decreasing order, thus breaking the spatial dependence. We then compared the bins of the resulting histogram in the data and in the simulated path.Mobile transectsFirst, the movement of the animal was linearly interpolated from the GPS data, meaning that between two recorded locations the individual followed a linear path. The speed of the animal between two locations was accordingly reconstructed using the recorded times (t_i) between each location. Second, we used mobile transects as the ecological sampling method, where each transect ‘count’ the intersection between its path and the animal’s one. The mobile transects followed a predefined path at a given constant speed as time increased. The area of vision of each transect was defined as a circle of a given radius. Each time the path of an individual collided with an area of vision, the count of the corresponding transect increased by 1. Two types of movements were used: linear and clockwise rotational transects. The initial locations of both types of transects are (X_1) and (X_F). Both the animal and mobile transects started to move at the same time. At each of the two locations (X_1), (X_F), 8 linear transects moved in the 8 cardinal directions, totalizing 16 transects. For the linear transects, every 10,000 time steps, we set (2 times 8) new transects starting at the same locations and following the same directions. Clockwise rotational transects were rotated around (X_1) and (X_F) using a 500 m radius. When we reached (t_n), we gathered the total count (i.e. the count of all transects). For the two types of transects, we gathered the total count for 6 different lines of sight: 50, 100, 200, 400, 500, 1000 m. and 4 speeds: s/4, s/2, s, 2 s with s the average speed of the animal. We then aggregated the overall count in each of the two types of transects, and compared the results from the data and the simulated path (Supplementary Methods and Fig. S2).Scale invarianceSeveral authors pointed out that the temporal resolution of the discretization is of importance: it should be relevant to the considered behavioral mechanisms5,48,49,50. Schlägel and Lewis focused on the quantification of movement models’ robustness under subsampled movement paths49. They found that increased subsampling leads to a strong deviation of the central parameter in resource selection models49,51. They underlined that important quantities derived from empirical data (e.g. parameters estimates, travel distance or sinuosity) can differ based on the temporal resolution of the data49,51. Moreover, Postlethwaite and Dennis highlighted the difficulty of comparing model results amongst tracking-datasets that vary substantially in temporal grain50). Each of the studied dataset has a relatively high sampling rate (roughly 10 m) and a period of study that is appropriate to the analysis of animal movement at the year scale (Table 1). In order to investigate such a possible effect on the BCR dynamic, we changed the sampling rate of the movement path to ensure that the three parameters (p_I), (p_s) and (p_F) are scale invariant. The movement path formed by the GPS observations (X_i) was subsampled (decimated) for each individual. We only kept every (k^{text {th}}) observation starting with the first one and (k in left[ 1,10right] ). For (k=1) the path corresponded to the original one. The time spent between each successive observation was also accordingly reconstructed in order to keep track of ({overline{T}}) in subsampled movement paths. The time between two locations (X_i) and (X_{i+k}) was reconstructed as:$$begin{aligned} t_{j}’ = sum _{i=j}^{i+(k-1)} t_i end{aligned}$$
(19)
with (j in left[ 1, 1+k, 1+2k, ldots , n-left( k-1 right) right] ). We did not change the value of (d_{min }) as we subsampled the movement path because we designed (p_s) for capturing GPS noise and movements that are associated with peculiar ecological behaviors that are beyond the scope of this study in terms of time and spatial scales (foraging for instance). We then compared the resulting parameters (p_I), (p_F) and (p_s) as the resampling rate k increased.FluctuationsWhereas the BCR is a stochastic process, the deterministic aspects of the 5 statistics were tested with an increasing number of steps (n_s). The statistic associated with each realization of the model (a simulated path) is a random variable. If the distribution of these random variables has low concentration (high variance) then it is not a convenient statistic as it cannot be used as a reference for assessing the model’s performance, even when averaging over multiples realizations. On the opposite, if the statistic is deterministic (no fluctuations) it can provide a reliable tool to assess the model’s performance. This was numerically tested over a range of increasing (n_s) values with (n_s = 10^4, 2times 10^4, ldots , 4 times 10^5). For each of those step values, a set of 100 BCR was simulated with parameters (p_I), (p_F) and (p_s) estimated from the first deer (see Table 2) and we studied the variance of the statistics.Sensitivity analysisIn order to assess whether the estimated parameters are optimal (i.e. providing the best possible performance) and to study parameter scarcity, we also evaluated the performance of the model using arbitrary weight values. We first started by evaluating how using one parameter instead of the three could alter the performance of the model. We then extended this sensitivity analyse by drawing arbitrary values for each parameter from a multi-dimensional square mesh, whose center corresponds to the estimated values of (p_I), (p_s), (p_F), estimated using GPS data (Fig. 1). We additionally used values that are distant from the estimated ones, up to (p_I=3), (p_F=3) and (p_s=5). We tested a total of 151 new configurations with these arbitrary values. For each configuration, we ran 150 simulations and evaluated them using the 5 statistics. The mean error of (|smash {{mathscr {S}}}- {tilde{mathscr {S}}}_k|) and its standard deviation are gathered and plotted for each arbitrary configuration. As a resume, we replicate the framework described in Fig. 1 but we inject arbitrary parameters instead of using data-driven parameterisations.ApplicationThe proposed model could be used to infer environmental and behavior information from the dataset. We chose to illustrate such an application by trying to detect anomalous voids (or holes) in the spatial territory of the individual using the GPS dataset and Monte-Carlo simulations of the model. Anomalous means that the observed void is not related to the randomness of the movement, but rather related to a geographical artifact. The parameters (p_I), (p_F), (p_s), (mu ) and (sigma ^2) of the BCR were accordingly estimated from the data of each individual, similarly to previous experiments (Fig. 1). A simple heuristic was used to find voids in empirical and simulated paths for each individual: we computed the alpha shape of all locations using a fixed alpha radius of 60 m. This allowed for determining the surface covered by all locations while preserving the voids. We then collected the area of each void provided they had an area of at least 100 m2. We focused on voids near the center of the alpha shape in order to avoid artificial voids, generated by the weak density of locations at the boundaries. We ran 10,000 iterations of the model for each animal and estimated the probability (p_{varnothing }) of finding voids of different sizes in the simulated paths. This probability was then compared to voids found in the GPS datasets and available environmental information was used to determine whether any geographical element(s) could explain the unexpected voids. More

1.Cheruvelil, K. S. & Soranno, P. A. Data-intensive ecological research is catalyzed by open science and team science. Bioscience 68, 813–822 (2018).Article 

Google Scholar 
2.Kelling, S. et al. Data-intensive science: a new paradigm for biodiversity studies. Bioscience 59, 613–620 (2009).Article 

Google Scholar 
3.Queiroz, N. et al. Global spatial risk assessment of sharks under the footprint of fisheries. Nature 572, 461–466 (2019).ADS 
CAS 
Article 

Google Scholar 
4.Clarke, R. Big data, big risks. Inf. Syst. J. 26, 77–90 (2016).Article 

Google Scholar 
5.Kroodsma, D. A. et al. Tracking the global footprint of fisheries. Science 359, 904–908 (2018).ADS 
CAS 
Article 

Google Scholar 
6.Braccini, M., Molony, B. & Blay, N. Patterns in abundance and size of sharks in northwestern Australia: cause for optimism. ICES J. Mar. Sci. 77, 72–82 (2020).Article 

Google Scholar 
7.Larcombe, J., Caton, A., Wiliams, D. M. & Speare, P. Western Tuna and Billfish Fisheries Research (Bureau of Resource Sciences, Canberra 1997).8.Larcombe, J. & Begg, G. Fishery Status Reports 2007. Status of Fish Stocks Managed by the Australian Government (Bureau of Rural Sciences, Canberra 2007).9.Hobsbawn, P. I., Patterson, H. M. & Blake, S. A. P. Australian National Report to the Scientific Committee of the Indian Ocean Tuna Commission (ABARES, Canberra 2020).10.Kaplan, R. M., Chambers, D. A. & Glasgow, R. E. Big data and large sample size: a cautionary note on the potential for bias. Clin. Transl. Sci. 7, 342–346 (2014).Article 

Google Scholar 
11.Bayraktarov, E. et al. Do big unstructured biodiversity data mean more knowledge? Front. Ecol. Evol. 6, 239 (2019).Article 

Google Scholar 
12.Braccini, M., Blay, N., Hesp, A. & Molony, B. Resource Assessment Report Temperate Demersal Elasmobranch Resource of Western Australia. Fisheries Research Report No. 294 (Department of Primary Industries and Regional Development, Perth, 2018).13.Roff, G., Brown, C. J., Priest, M. A. & Mumby, P. J. Decline of coastal apex shark populations over the past half century. Commun. Biol. 1, 223 (2018).Article 

Google Scholar 
14.Hampton, S. E. et al. Big data and the future of ecology. Front. Ecol. Environ. 11, 156–162 (2013).Article 

Google Scholar 
15.R Core Team. R: A language and environment for statistical computing. http://www.R-project.org/ (R Foundation for Statistical Computing, 2020). More