Ecology

Lehtonen12 presents three simple models with the same broad structure: a single mutant individual with divergent mating behaviour arises in a population of ‘residents’ that all play the same strategy, and the success of that mutant is then followed (Figs. 1, 2). Specifically, Lehtonen investigates the fitness benefits of increased mating for mutant males in comparison to mutant females. Two important parameters can be varied: (i) the degree of anisogamy (defined here as the ratio of sperm number to egg number), which captures how divergent males and females are in the size (and thus number) of gametes they produce, and (ii) the efficiency of fertilisation, which determines how easily gametes can find and fuse with each other. If fertilisation is highly efficient, then gametes of the less numerous type will achieve nearly full fertilisation; on the other hand, inefficient fertilisation can result in gametes of both sexes going unfertilised.Fig. 2: Structure of the three models of Lehtonen12, showing differences in mating behaviour between resident males (green), resident females (blue) and mutant males and females (both yellow).For illustration, we suppose that females produce four eggs each and males produce eight sperm (the anisogamy ratio in nature is typically much higher). In Model 1, resident individuals spawn monogamously in a ‘nest’ (black outline), whereas mutant males and females can bring additional partners to their nest to spawn in a group. In Model 2, resident individuals divide their gametes equally among m spawning groups, each consisting of m individuals of each sex (shown here with m = 2). Mutant males and females instead divide their gametes among a larger or smaller number of groups, mmutant (shown here with mmutant = 4). In Model 3, there is a further sex asymmetry in addition to anisogamy: Fertilisation takes place inside the female’s body. Resident individuals mate with m partners (shown here with m = 2), whereas mutant males and females mate with a larger or smaller number of partners, mmutant (shown here with mmutant = 4).Full size imageIn the first two models, fertilisation is external and no assumptions are made about pre-existing differences between the sexes apart from the number of gametes they produce. In other words, males and females are identical except that males produce sperm in greater numbers than females produce eggs. In Model 1, resident individuals are assumed to mate monogamously, whereas a mutant can monopolise multiple partners of the opposite sex (Fig. 2). Importantly, both male and female mutants can bring additional partners back to their ‘nest’ to spawn in a group. When fertilisation is highly efficient, females can fertilise all of their eggs by bringing back a single male, and there is simply no benefit (in this model) of seeking further partners (Fig. 1A). In contrast, anisogamy means that males always produce at least some gametes in excess, and thus can benefit from seeking additional mates. When fertilisation is inefficient, however, both sexes benefit from increasing the concentration of opposite-sex gametes at their ‘nest’ (Fig. 1B). This latter benefit is sex-symmetric, whereas the former continues to apply only to males. As a consequence, the Bateman gradients are always steeper for males than for females (Fig. 1A, B), confirming Bateman’s argument.Model 2 similarly assumes external fertilisation, but in this case the resident males and females meet in groups consisting of m individuals of each sex (Fig. 2). Fertilisation occurs via group spawning. It is assumed that each resident individual divides its gametes evenly across M groups, whereas mutant individuals can instead spread their gametes over a larger or smaller number of groups (note that the author assumes that M = m, but this assumption could be relaxed without undermining the core argument). Spreading gametes out across a larger number of spawning groups does not increase the concentration of opposite-sex gametes they encounter (Fig. 2). However, a mutant that spreads its gametes more widely reduces the density of its own gametes across those groups in which it spawns. This in turn results in there being more opposite-sex gametes for each gamete of the mutant’s sex in those groups. For example, in Fig. 2, mutant males spawn in twice as many groups as resident males and thereby halve the density of their own sperm in each group. The resulting egg-to-sperm ratio of (frac{4}{6}=frac{2}{3}) is more favourable than the ratio of (frac{4}{8}=frac{1}{2}) that the resident males experience. Mutant females can similarly increase local sperm-to-egg ratios by spreading their eggs over more groups. However, in contrast to males, this only leads to fitness benefit if fertilisation is inefficient, and even then the benefit to females is very modest (scarcely perceptible in Fig. 1D). Gamete spreading reduces wasteful competition among the mutants’ own gametes for fertilisation. Such ‘local’ gamete competition, like gamete competition more generally, is stronger among sperm than among eggs because sperm are more numerous under anisogamy13,14. Consequently, as in Model 1, Bateman gradients are always steeper in males (Fig. 1C, D). Recall that the results of the above models emerge in the absence of any assumptions beyond the sex difference in the number of gametes produced.The third and final model allows for a further pre-existing difference between the sexes in addition to anisogamy: internal fertilisation, which is common and widespread in animals (Fig. 2)15. Each female is assumed to mate with m males, while each male divides his gametes evenly among m females. As in the previous two models, males benefit more than females from additional matings under most conditions. However, in the particular case where fertilisation is highly inefficient and the ratio of sperm to eggs is not too large, the pattern can theoretically reverse, such that female Bateman gradients exceed their male counterparts (Fig. 1F). The reason is that the effects of gamete concentration are asymmetric under internal fertilisation: Multiple mating by a female increases the local concentration of sperm its eggs experience, whereas a male’s multiple mating does not increase the concentration of eggs around its sperm (Fig. 2). Under conditions of severe sperm limitation—due to both weak anisogamy and highly inefficient fertilisation—this can lead to females benefitting more from additional matings than males (Fig. 1F). Although intriguing, it is unclear whether this finding has any empirical relevance, as sperm limitation is probably rarely severe in internal fertilisers. Under more realistic conditions of moderate to high fertilisation rates, sex differences in the degree of local gamete competition once again become decisive, and male Bateman gradients exceed their female counterparts (Fig. 1E). More

Hauer, M. E. et al. Sea-level rise and human migration. Nat. Rev. Earth Environ. 1, 28–39 (2020).ADS 
Article 

Google Scholar 
Duy, P., Chapman, L., Tight, M., Thuong, L. & Linh, P. Urban resilience to floods in coastal cities: Challenges and opportunities for Ho Chi Minh city and other emerging cities in southeast Asia. J. Urban Plan. Dev. 144, 05017018 (2018).Article 

Google Scholar 
Magno, R. et al. Semi-automatic operational service for drought monitoring and forecasting in the Tuscany region. Geosciences 8, 49 (2018).ADS 
Article 

Google Scholar 
Rico, A., Olcina, J., Baños, C., Garcia, X. & Sauri, D. Declining water consumption in the hotel industry of mass tourism resorts: Contrasting evidence for Benidorm, Spain. Curr. Issues Tour. 23, 770–783 (2020).Article 

Google Scholar 
Hasnat, G. T., Kabir, M. A. & Hossain, M. A. Major environmental issues and problems of South Asia, particularly Bangladesh. Handb. Environ. Mater. Manag., 1–40 (2018).Neumann, B., Vafeidis, A. T., Zimmermann, J. & Nicholls, R. J. Future coastal population growth and exposure to sea-level rise and coastal flooding—A global assessment. PLoS One 10, e0118571 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Cao, A. et al. Future of Asian Deltaic Megacities under sea level rise and land subsidence: Current adaptation pathways for Tokyo, Jakarta, Manila, and Ho Chi Minh City. Curr. Opin. Environ. Sustain. 50, 87–97 (2021).Article 

Google Scholar 
Rahmasary, A. N. et al. Overcoming the challenges of water, waste and climate change in Asian cities. Environ. Manag. 63, 520–535 (2019).ADS 
Article 

Google Scholar 
Smol, M., Adam, C. & Preisner, M. Circular economy model framework in the European water and wastewater sector. J. Mater. Cycles Waste Manag. 22, 682–697 (2020).Article 

Google Scholar 
Islam, M. F., Bhattacharya, B. & Popescu, I. Flood risk assessment due to cyclone-induced dike breaching in coastal areas of Bangladesh. Nat. Hazards Earth Syst. Sci. 19, 353–368 (2019).ADS 
Article 

Google Scholar 
Salim, M. A. & Siswanto, A. B. Kajian Penanganan Dampak Banjir Kabupaten Pekalongan. Rang Tek. J. 4, 295–303 (2021).Article 

Google Scholar 
Endo, A. et al. Describing and visualizing a water–energy–food nexus system. Water 10, 1245 (2018).Article 

Google Scholar 
Gurdak, J. J., Geyer, G. E., Nanus, L., Taniguchi, M. & Corona, C. R. Scale dependence of controls on groundwater vulnerability in the water–energy–food nexus, California Coastal Basin aquifer system. J. Hydrol. Reg. Stud. 11, 126–138 (2017).Article 

Google Scholar 
Lu, J., Lin, Y., Wu, J. & Zhang, C. Continental-scale spatial distribution, sources, and health risks of heavy metals in seafood: Challenge for the water-food-energy nexus sustainability in coastal regions?. Environ. Sci. Pollut. Res. 28, 63815–63828 (2021).CAS 
Article 

Google Scholar 
Miller-Robbie, L., Ramaswami, A. & Amerasinghe, P. Wastewater treatment and reuse in urban agriculture: Exploring the food, energy, water, and health nexus in Hyderabad, India. Environ. Res. Lett. 12, 075005 (2017).ADS 
Article 
CAS 

Google Scholar 
Taniguchi, M., Endo, A., Gurdak, J. J. & Swarzenski, P. Water-energy-food nexus in the Asia-Pacific region. J. Hydrol. 11, 1–8 (2017).
Google Scholar 
Bahri, M. Analysis of the water, energy, food and land nexus using the system archetypes: A case study in the Jatiluhur reservoir, West Java, Indonesia. Sci. Total Environ. 716, 137025 (2020).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Lubis, R., Delinom, R., Martosuparno, S. & Bakti, H. Water-Food Nexus in Citarum Watershed, Indonesia Vol. 118, 012023 (IOP Publishing, 2018).
Google Scholar 
Pawitan, H., Delinom, R. & Taniguchi, M. The human–environment sustainability in Indonesia: The case of the Citarum basin Vol. 23 (UNESCO-IHP, 2015).Carmichael, L. et al. Urban planning as an enabler of urban health: Challenges and good practice in England following the 2012 planning and public health reforms. Land Use Policy 84, 154–162 (2019).Article 

Google Scholar 
World Health Organization. Addressing the Social Determinants of Health: The Urban Dimension and the Role of Local Government (World Health Organization, 2012).
Google Scholar 
Trencher, G. & Karvonen, A. Stretching, “smart”: Advancing health and well-being through the smart city agenda. Local Environ. 24, 610–627 (2019).Article 

Google Scholar 
Yang, L. et al. Can an island economy be more sustainable? A comparative study of Indonesia, Malaysia, and the Philippines. J. Clean. Prod. 242, 118572 (2020).Article 

Google Scholar 
Choirunisa, A. K. & Giyarsih, S. R. Kajian Kerentanan Fisik, Sosial, dan Ekonomi Pesisir Samas Kabupaten Bantul Terhadap Erosi Pantai. J. Bumi Indones. 5 (2016).Gumay, A. Validity and reliability maritime English seafarers proficiency test. INFERENCE J. Engl. Lang. Teach. 3, 64–69 (2021).Article 

Google Scholar 
Tarigan, M. S. Perubahan garis pantai di wilayah pesisir perairan Cisadane, Provinsi Banten. Makara J. Sci. (2010).Pruss-Ustun, A., Corvalán, C. F., World Health Organization. Preventing Disease Through Healthy Environments: Towards an Estimate of the Environmental Burden of Disease (World Health Organization, 2006).
Google Scholar 
Baasanjargal, T., Soon-Joo, A. & Mi-Jeong, K. Comparative analysis of Indonesian Batik traditional patterns: Focused on patterns of Yogyakarta and Pekalongan in Java Island. 한복문화 22, 75–91 (2019).Rismawati, S. D., Sofiani, T. & Rahmawati, D. R. Legal culture of religious capitalism on Batik business (a case study in Pekalongan Indonesia). JL Pol. Glob. 33, 107 (2015).
Google Scholar 
Pekalongan, B. K. Kota Pekalongan dalam Angka 2021 (2021).Google Maps. Pekalongan, Central Java (2022).Sunarjo, W. A., Ilmiani, A. & Ardianingsih, A. Analisis SWOT Sebagai Pengembangan UMKM Berbasis Ekonomi Kreatif Destinasi Pariwisata Batik Kota Pekalongan. Pena J. Ilmu Pengetah. Dan Teknol. 33, 34–43 (2019).Article 

Google Scholar 
Perpustakaan Provinsi Jawa Tengah. Museum Batik Pekalongan (2017).Brzezina, N. et al. Development of organic farming in Europe at the crossroads: Looking for the way forward through system archetypes lenses. Sustainability 9, 821 (2017).Article 

Google Scholar 
Gillies, A. & Maliapen, M. Using healthcare system archetypes to help hospitals become learning organisations. J. Model. Manag. (2008).Braun, W. The System Archetypes. The Systems Modeling Workbook, 1–26 (2002).Sterman, J. System Dynamics: Systems thinking and modeling for a complex world (2002).Islam, M. & Raja, D. R. Waterlogging risk assessment: An undervalued disaster risk in coastal urban community of Chattogram, Bangladesh. Earth 2, 151–173 (2021).Article 

Google Scholar 
Brzezina, N., Kopainsky, B. & Mathijs, E. Can organic farming reduce vulnerabilities and enhance the resilience of the European food system? A critical assessment using system dynamics structural thinking tools. Sustainability 8, 971 (2016).Article 

Google Scholar 
Nguyen, N. C. & Bosch, O. J. A systems thinking approach to identify leverage points for sustainability: A case study in the Cat Ba Biosphere Reserve, Vietnam. Syst. Res. Behav. Sci. 30, 104–115 (2013).Article 

Google Scholar 
Maani, K. E. & Cavana, R. Y. Systems Thinking, System Dynamics: Managing Change and Complexity (Pearson Prentice Hall, 2007).
Google Scholar 
Braun, W. The System Archetypes-the Systems Modeling Workbook. Available Wwwu Uniklu Ac Atgossimitpapsdwbsysarch Pdf (2002).Senge, P. M. The Fifth Discipline: The Art and Practice of the Learning Organization (Currency, 2006).
Google Scholar 
Bahri, M. et al. Deliverable 3.3: Integrated model with ad-hoc systems model of urban water supply (2018).Pekalongan, B. P. P. D. K. Pekalongan dalam Angka (2021).Fajar, M., Mediani, A. & Finesa, Y. Analisis Peranan IPAL dalam Strategi Penanganan Limbah Industri Batik di Kota Pekalongan. in Prosiding Seminar Nasional Geografi UMS X 2019 (2019).Kartika, F. D. S. & Helmi, M. Meta-analysis of Community’s Adaptation Pattern with Tidal Flood in Pekalongan City, Central Java, Indonesia Vol. 125, 09001 (EDP Sciences, 2019).
Google Scholar 
Kartika, F. D. S., Helmi, M. & Amirudin, A. Analisis Perubahan Penggunaan Lahan di Wilayah Pesisir Kota Pekalongan Menggunakan Citra Lansat 8, vol. 1 (2019).Damayanti, M. & Latifah, L. Strategi Kota Pekalongan dalam pengembangan wisata kreatif berbasis industri batik. J. Pengemb. Kota 3, 100–111 (2017).Article 

Google Scholar 
Pekalongan, B. K. Kota Pekalongan dalam Angka 2002 (2002).Andreas, H., Abidin, H. Z., Sarsito, D. A. & Pradipta, D. Adaptation of ‘early climate change disaster’ to the Northern coast of Java Island Indonesia. Eng. J. 22, 207–219 (2018).Article 

Google Scholar 
Marfai, M. A. et al. The impact of tidal flooding on a coastal community in Semarang, Indonesia. Environmentalist 28, 237–248 (2008).Article 

Google Scholar 
Chaussard, E., Amelung, F., Abidin, H. & Hong, S.-H. Sinking cities in Indonesia: ALOS PALSAR detects rapid subsidence due to groundwater and gas extraction. Remote Sens. Environ. 128, 150–161 (2013).ADS 
Article 

Google Scholar 
Andreas, H., Abidin, H. Z., Sarsito, D. A. & Pradipta, D. Remotes Sensing Capabilities on Land Subsidence and Coastal Water Hazard and Disaster Studies Vol. 500, 012036 (IOP Publishing, 2020).
Google Scholar 
Shofiana, R., Subardjo, P. & Pratikto, I. Analisis perubahan penggunaan lahan di wilayah pesisir Kota pekalongan menggunakan data landsat 7 etm+. J. Mar. Res. 2, 35–43 (2013).
Google Scholar 
Wijaya, A. Analisis Dinamika Pola Spasial Penggunaan Lahan Pada Wilayah Terdampak Kenaikan Muka Air Laut di Kota Pekalongan (2017).El-Fath, D. D. I., Atmodjo, W., Helmi, M., Widada, S. & Rochaddi, B. Analisis Spasial Area Genangan Banjir Rob Setelah Pembangunan Tanggul di Kabupaten Pekalongan, Jawa Tengah. Indones. J. Oceanogr. 4, 96–110 (2022).
Google Scholar 
Novita, M. G., Helmi, M., Widiaratih, R., Hariyadi, H. & Wirasatriya, A. Mengkaji Area Genangan Banjir Pasang Terhadap Penggunaan Lahan Pesisir Tahun 2020 Menggunakan Metode Geospasial di Kabupaten Pekalongan, Provinsi Ja. Indones. J. Oceanogr. 3, 14–26 (2021).Article 

Google Scholar 
Salim, M. A. Penanganan Banjir dan Rob di Wilayah Pekalongan. J. Tek. Sipil 11, 15–23 (2018).
Google Scholar 
Jumatiningrum, N. & Indrayati, A. Strategi Adaptasi Masyarakat Kelurahan Bandengan Kecamatan Pekalongan Utara dalam Menghadapi Banjir Pasang Air Laut (Rob). Edu Geogr. 9, 136–143 (2021).
Google Scholar 
BNPB. Data Kebencanaan Nasional (BNPB, 2021).
Google Scholar 
Giampietro, M., Aspinall, R. J., Ramos-Martin, J. & Bukkens, S. G. Resource Accounting for Sustainability Assessment: The Nexus Between Energy, Food, Water and Land Use (Routledge, 2014).Book 

Google Scholar 
Meadows, D. H., Randers, J. & Meadows, D. L. The Limits to Growth (1972) (Yale University Press, 2013).MATH 

Google Scholar 
Albrecht, T., Crootof, A. & Scott, C. Trends in the development of water–energy–food nexus methods (2017).Leck, H., Fitzpatrick, D. & Burchell, K. Energy, water and food: Towards a critical nexus approach. in Handbook on the Geographies of Energy (Edward Elgar Publishing, 2017).Scott, C. A., Kurian, M. & Wescoat, J. L. The water–energy–food nexus: Enhancing adaptive capacity to complex global challenges. in Governing the Nexus 15–38 (Springer, 2015).Wanty, E. E. Analisis Produksi Batik Cap Dari UKM Batik Kota Pekalongan (Studi Pada Sentra Batik Kota Pekalongan-Jawa Tengah, 2006).
Google Scholar 
Mankiw, N. G. Macroeconomics Vol. 41 (Worth Publishers, 2003).
Google Scholar 
Shen, J. & Kee, G. Development and Planning in Seven Major Coastal Cities in Southern and Eastern China (Springer, 2017).Book 

Google Scholar 
Xu, C., Haase, D., Su, M. & Yang, Z. The impact of urban compactness on energy-related greenhouse gas emissions across EU member states: Population density vs physical compactness. Appl. Energy 254, 113671 (2019).Article 

Google Scholar 
Marfai, M. A. & Cahyadi, A. Dampak bencana banjir pesisir dan adaptasi masyarakat terhadapnya di kabupaten Pekalongan (2017).Wartadesa.net. Tiga hari banjir rendam Pekalongan (2018).Google Maps. A dike in Pekalongan (n.d).Anindita, R. M., Susilowati, I. & Muhammad, F. Analisis Efektifitas Tanggul Laut di Pesisir Pekalongan Terhadap Penurunan Intensitas Banjir, vol. 2 80–88 (2020).Taniguchi, M. Groundwater and Subsurface Environments: Human Impacts in Asian Coastal Cities (Springer Science & Business Media, 2011).Book 

Google Scholar 
Baños, C. J., Hernández, M., Rico, A. M. & Olcina, J. The hydrosocial cycle in coastal tourist destinations in Alicante, Spain: Increasing resilience to drought. Sustainability 11, 4494 (2019).Article 

Google Scholar 
Sauda, R. H. & Nugraha, A. L. Kajian pemetaan kerentanan banjir rob di kabupaten pekalongan. J. Geod. Undip 8, 466–474 (2019).
Google Scholar 
Wartadesa.net. Ratusan warga Sragi masih mengungsi (2022).Buchori, I. et al. Adaptation to coastal flooding and inundation: Mitigations and migration pattern in Semarang City, Indonesia. Ocean Coast. Manag. 163, 445–455 (2018).Article 

Google Scholar 
Setiadi, R. & Nalau, J. Can urban regeneration improve health resilience in a changing climate? (2015).Isham, A., Mair, S. & Jackson, T. Wellbeing and productivity: A review of the literature (2020).Banson, K. E., Nguyen, N. C. & Bosch, O. J. Using system archetypes to identify drivers and barriers for sustainable agriculture in Africa: A case study in Ghana. Syst. Res. Behav. Sci. 33, 79–99 (2016).Article 

Google Scholar 
Lavrnić, S., Zapater-Pereyra, M. & Mancini, M. Water scarcity and wastewater reuse standards in Southern Europe: Focus on agriculture. Water. Air Soil Pollut. 228, 1–12 (2017).Article 
CAS 

Google Scholar 
Tortajada, C. & Nam Ong, C. Reused water policies for potable use (2016).Murali, R. M., Riyas, M., Reshma, K. & Kumar, S. S. Climate change impact and vulnerability assessment of Mumbai city, India. Nat. Hazards 102, 575–589 (2020).Article 

Google Scholar 
Abdullah, A. Y. M. et al. Spatio-temporal patterns of land use/land cover change in the heterogeneous coastal region of Bangladesh between 1990 and 2017. Remote Sens. 11, 790 (2019).ADS 
Article 

Google Scholar 
Ginanjar, A., Rezagama, A. & Handayani, D. S. Rencana Induk Sistem Penyediaan Air Minum Kota Pekalongan (2015).Reiblich, J., Hartge, E., Wedding, L., Killian, S. & Verutes, G. Bridging climate science, law, and policy to advance coastal adaptation planning. Mar. Policy 104, 125–134 (2019).Article 

Google Scholar 
Cook, B. I. et al. Revisiting the leading drivers of Pacific coastal drought variability in the contiguous United States. J. Clim. 31, 25–43 (2018).ADS 
Article 

Google Scholar 
Jodar-Abellan, A., Valdes-Abellan, J., Pla, C. & Gomariz-Castillo, F. Impact of land use changes on flash flood prediction using a sub-daily SWAT model in five Mediterranean ungauged watersheds (SE Spain). Sci. Total Environ. 657, 1578–1591 (2019).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Thanvisitthpon, N., Shrestha, S. & Pal, I. Urban flooding and climate change: A case study of Bangkok, Thailand. Environ. Urban. Asia 9, 86–100 (2018).Article 

Google Scholar 
Laksmi, G. S. Dampak Alih Fungsi Lahan dan Curah Hujan terhadap Banjir di Kota Pekalongan, Jawa Tengah, 382–391 (2020).Dhiman, R., VishnuRadhan, R., Eldho, T. & Inamdar, A. Flood risk and adaptation in Indian coastal cities: Recent scenarios. Appl. Water Sci. 9, 1–16 (2019).ADS 
Article 

Google Scholar 
Bahri, M. & Cremades, R. The Urban Drought Nexus Tool. Zenodo (2021). More

Eco-evolutionary model on spatial graphsWe establish an individual-based model (IBM) where individuals are structured over a trait space and a graph representing a landscape. For the sake of simplicity, we consider the case of asexual reproduction and haploid genetics29. Individuals die, reproduce, mutate and migrate in a stochastic fashion, which together results in macroscopic properties. The formulation of the stochastic IBM allows an analytical description of the dynamics at the population level, which links emergent properties to the elementary processes that generate them.The trait space ({{{{{{{mathcal{X}}}}}}}}subseteq {{mathbb{R}}}^{d}) is continuous and can be split into a neutral trait space ({{{{{{{mathcal{U}}}}}}}}) and an adaptive trait space ({{{{{{{mathcal{S}}}}}}}}). We refer to neutral traits (uin {{{{{{{mathcal{U}}}}}}}}) as traits that are not under selection, in contrast to adaptive traits (sin {{{{{{{mathcal{S}}}}}}}}), which experience selection. The graph denoted by G is composed of a set of vertices {v1,v2,…,vM} that correspond to habitat patches (suitable geographical areas), and a set of edges that constrain the movement of individuals between the habitat patches. We use the original measure of genetic differentiation for quantitative traits QST (standing for Q-statistics) in the case of haploid populations45,46. We denote the neutral trait value of the kth individual on vi as ({u}_{k}^{(i)}), the number of individuals on vi as N(i), the mean neutral trait on vi as ({overline{u}}^{(i)}), and the mean neutral trait in the metapopulation as (overline{u}). It follows that we quantify neutral differentiation QST,u as$${Q}_{ST,u}={sigma }_{B,u}^{2}/({sigma }_{B,u}^{2}+{sigma }_{W,u}^{2})$$
(1)
where ({sigma }_{B,u}^{2}={mathbb{E}}[frac{1}{M}{sum }_{i}{left({overline{u}}^{(i)}-overline{u}right)}^{2}]) denotes the expected neutral trait variance between the vertices and ({sigma }_{W,u}^{2}=frac{1}{M}mathop{sum }nolimits_{i}^{M}{mathbb{E}}left[frac{1}{{N}^{(i)}}{sum }_{k}{left({u}_{k}^{(i)}-{overline{u}}^{(i)}right)}^{2}right]) denotes the average expected neutral trait variance within vertices. We similarly quantify adaptive differentiation QST,s.Following the Gillespie update rule47, individuals with trait ({x}_{k}in {{{{{{{mathcal{X}}}}}}}}) on vertex vi are randomly selected to give birth at rate b(i)(xk) and die at rate d(N(i)) = N(i)/K, where K is the local carrying capacity. The definition of d therefore captures competition, which is proportional to the number of individuals on a vertex and does not depend on the individuals’ traits (we relax this assumption later on). The offspring resulting from a birth event inherits the parental traits, which can independently be affected by mutations with probability μ. A mutated trait differs from the parental trait by a random change that follows a normal distribution with variance ({sigma }_{mu }^{2}) (corresponding to the continuum of alleles model48). The offspring can further migrate to neighbouring vertices by executing a simple random walk on G with probability m. A schematic overview of the two different settings considered is provided in Fig. 1. Under the setting with no selection, individuals are only characterised by neutral traits so that ({{{{{{{mathcal{X}}}}}}}}={{{{{{{mathcal{U}}}}}}}}). For individuals on a vertex with trait xk ≡ uk we define b(i)(xk) ≡ b, so that the birth rate is constant. This ensures that neutral traits do not provide any selective advantage. Under the setting with heterogeneous selection, each vertex of the graph vi is labelled by a habitat type with environmental condition Θi that specifies the optimal adaptive trait value on vi. It follows that, for individuals with traits ({x}_{k}=({u}_{k},{s}_{k})in {{{{{{{mathcal{U}}}}}}}}times {{{{{{{mathcal{S}}}}}}}}) on vi, we define$${b}^{(i)}({x}_{k})equiv {b}^{(i)}({s}_{k})=b(1-p{({s}_{k}-{{{Theta }}}_{i})}^{2})$$
(2)
where p is the selection strength41. This ensures that the maximum birth rate on vi is attained for sk = Θi, which results in a differential advantage that acts as an evolutionary stabilising force. In the following we consider two habitat types denoted by I and II with symmetric environmental conditions θI and θII, so that Θi ∈ {θI, θII} and θII = − θI = θ, where θ can be viewed as the habitat heterogeneity41.Fig. 1: Graphical representation of the structure of individuals in the eco-evolutionary model.a Setting with no selection, where individuals are characterised by a set of neutral traits (uin {{{{{{{mathcal{U}}}}}}}}). The scatter plots represent a projection of the first two components of u for the individuals present on the designated vertices at time t = 1000, obtained from one simulation of the IBM. b Setting with heterogeneous selection. In this setting, individuals are additionally characterised by adaptive traits (sin {{{{{{{mathcal{S}}}}}}}}). Blue vertices favour the optimal adaptive trait value θI, while red vertices favour θII. The scatter plots represent a projection of the first component of u and s for the individuals present on the designated vertices at time t = 1000, obtained from one simulation. The majority of individuals are locally well-adapted and have an adaptive trait close to the optimal value, but some maladaptive individuals originating from neighbouring vertices are also present. m = 0.05.Full size imageDeterministic approximation of the population dynamics under no selectionThe model can be formulated as a measure-valued point process (30 and Supplementary Note). Under this formalism, we demonstrate in the Supplementary Note how the population size and the trait dynamics show a deterministic behaviour when a stabilising force dampens the stochastic fluctuations. This makes it possible to express the dynamics of the macroscopic properties with deterministic differential equations, connecting emergent patterns to the processes that generate them. In particular, in the setting of no selection, competition stabilises the population size fluctuations, and the dynamics can be considered deterministic and expressed as$${partial }_{t}{N}_{t}^{(i)}={N}_{t}^{(i)}left[b(1-m)-frac{{N}_{t}^{(i)}}{K}right]+mbmathop{sum}limits_{jne i}frac{{a}_{i,j}}{{d}_{j}}{N}_{t}^{(j)}$$
(3)
where (A={({a}_{i,j})}_{1le i,jle M}) is the adjacency matrix of the graph G and D = (d1,d2,…,dM) is a vector containing the degree of each vertex (number of edges incident to the vertex). The first term on the right-hand side corresponds to logistic growth, which accounts for birth and death events of non-migrating individuals. The second term captures the gains due to migrations, which depend on the graph topology. Assuming that all vertices with the same degree have an equivalent position on the graph, corresponding to a “mean field” approach (see Methods), one can obtain a closed-form solution from Eq. (3) (see Eq. (12)), which shows that the average population size (overline{N}) scales with ({langle sqrt{k}rangle }^{2}/langle krangle), where 〈k〉 is the average vertex degree and (langle sqrt{k}rangle) is the average square-rooted vertex degree. The quantity ({langle sqrt{k}rangle }^{2}/langle krangle), denoted as hd, relates to the homogeneity in vertex degree of the graph and can therefore be viewed as a measure negatively associated with heterogeneity in connectivity. Simulations of the IBM illustrate that hd can explain differences in population size for complex graph topologies with varying migration regimes (Fig. 2a for graphs with M = 7 vertices and Supplementary Fig. 1a for M = 9). This analytical result is connected to theoretical work on reaction-diffusion processes49 and highlights that irregular graphs (graphs whose vertices do not have the same degree) result in unbalanced migration fluxes that affect the ecological balance between births and deaths. Highly connected vertices present an oversaturated carrying capacity (N(i)  > bK, see Methods), increasing local competition and lowering total population size compared with regular graphs (Fig. 2a). Because populations with small sizes experience more drift (31 and Supplementary Fig. 2), this result indicates that graph topology affects neutral differentiation not only through population isolation, but also by affecting population dynamics.Fig. 2: Effect of and hd on average population size (overline{N}) and neutral differentiation QST,u in the setting with no selection.a Response of (overline{N}) to homogeneity in degree ({h}_{d}={langle sqrt{k}rangle }^{2}/langle krangle) for all undirected connected graphs with M = 7 vertices and m = 0.5. b Response of QST,u to average path length for similar simulations obtained with m = 0.01. c Response of QST,u to homogeneity in degree hd for the same data. In a, b, and c, each dot represents average results from 5 replicate simulations of the IBM, the colour scale corresponds to the proportion of the graphs with similar x and y-axis values (graph density), and the blue line corresponds to a linear fit. d Standardized effect of hd and on QST,u, obtained from multivariate regression models independently fitted on similar data obtained for m = 0.01 and m = 0.5. The contributions of and hd to QST,u are alike for low migration regimes. Error bars show 95% confidence intervals. Analogous results on graphs with M = 9 vertices are presented in Supplementary Fig. 1 and all regression details can be found in Supplementary Table 2.Full size imageNonetheless, the stochasticity of the processes at the individual level can propagate to the population level and substantially affect the macroscopic properties. In particular, neutral differentiation emerges from the stochastic fluctuations of the populations’ neutral trait distribution. These fluctuations complicate an analytical underpinning of the dynamics, and in this case simulations of IBM offer a straightforward approach to evaluate the level of neutral differentiation.Effect of graph topology on neutral differentiation under no selectionWe study a setting with no selection and investigate the effect of the graph topology on neutral differentiation. When migration is limited, individuals’ traits are coherent on each vertex but stochastic drift at the population level generates neutral differentiation between the vertices. Migration attenuates neutral differentiation because it has a correlative effect on local trait distributions. Following21,22,26, we expect that the intensity of the correlative effect depends on the average path length of the graph 〈l〉, defined as the average shortest path between all pairs of vertices50. For a constant number of vertices, 〈l〉 is strictly related to the mean betweenness centrality and quantifies the graph connectivity50. High 〈l〉 implies low connectivity and greater isolation of populations, and hence we expect that graphs with high 〈l〉 are associated with high differentiation levels. We consider various graphs with an identical number of vertices and run simulations of the IBM to obtain the neutral differentiation level QST,u attained after a time long enough to discard transient dynamics (see Methods). We then interpret the discrepancies in QST,u across the simulations by relating them to the underlying graph topologies.We observe strong differences in QST,u across graphs for varying m, and find that 〈l〉 explains at least 55% of the variation in QST,u across all graphs with M = 7 vertices for (Fig. 2b). Nonetheless, some specific graphs, such as the star graph, present higher levels of QST,u than expected by their average path length. To explain this discrepancy, we explore the effect of homogeneity in vertex degree hd, as we showed in Eq. (12) that it decreases population size, which should in turn increase QST,u by intensifying stochastic drift. We find that hd explains 57% of the variation for low m (Fig. 2c). However, the fit remains similar after correcting for differences in population size (see Supplementary Table 1), indicating that irregular graphs structurally amplify the isolation of populations. Unbalanced migration fluxes lead central vertices to host more individuals than allowed by their carrying capacity. This causes increased competition that results in a higher death rate, so that migrants have a lower probability of further spreading their trait. Highly connected vertices therefore behave as bottlenecks, increasing the isolation of peripheral vertices and consequently amplifying QST,u.We then evaluate the concurrent effect of 〈l〉 and hd on QST,u with a multivariate regression model that we fit independently for low and high migration regimes (Fig. 2d). The multivariate regression model explains at least 70% of the variation in QST,u for the migration regimes considered and for graphs with M = 7 vertices (see Supplementary Table 2 for details). Moreover, we find that 〈l〉 and hd have akin contributions to neutral differentiation for low m, but the effect of 〈l〉 increases for higher migration regimes while the effect of hd decreases. To ensure that these conclusions can be generalised to larger graphs, we conduct the same analysis on a subset of graphs with M = 9 vertices and find congruent results (Supplementary Fig. 1). In the absence of selection and with competitive interactions, graphs with a high average path length 〈l〉 and low homogeneity in vertex degree hd, or similarly graphs with low connectivity and high heterogeneity in connectivity, show high levels of neutral differentiation.Deterministic approximation of the population dynamics and adaptation under heterogeneous selectionWe next consider heterogeneous selection and investigate the response of adaptive differentiation to the spatial distribution of habitat types, denoted as the Θ-spatial distribution. Adaptive differentiation emerges from local adaptation, but migration destabilises adaptation as a result of the influx of maladaptive migrants. We expect that higher connectivity between vertices of similar habitat type increases the level of adaptive differentiation, because it increases the proportion of well-adapted migrants. Local adaptation can be investigated by approximating the stochastic dynamics of the trait distribution with a deterministic partial differential equation (PDE). We demonstrate under mean-field assumption how the deterministic approximation can be reduced to an equivalent two-habitat model. We analyse the reduced model with the theory of adaptive dynamics36,41 and find a critical migration threshold m⋆ that determines local adaptation. m⋆ depends on a quantity coined the habitat assortativity rΘ, and we demonstrate with numerical simulations that rΘ determines the overall adaptive differentiation level QST,s reached at steady state in the deterministic approximation.Heterogeneous selection, captured by the dependence of the birth rate on Θi, generates a stabilising force that dampens the stochastic fluctuations of the adaptive trait distribution. The dynamics of the adaptive trait distribution consequently shows a deterministic behavior and we demonstrate in the Supplementary Note and Supplementary Figs. 3 and 4 that the number of individuals on vi with traits (sin {{Omega }}subset {{{{{{{mathcal{S}}}}}}}}) can be approximated by the quantity ∫Ωn(i)(s)ds, where n(i) is a continuous function solution of the PDE$${partial }_{t}{n}_{t}^{(i)}(s)= , {n}_{t}^{(i)}(s)left[{b}^{(i)}(s)(1-m)-frac{1}{K}{int}_{{{{{{{{mathcal{S}}}}}}}}}{n}_{t}^{(i)}({{{{{{{bf{s}}}}}}}})d{{{{{{{bf{s}}}}}}}}right]\ +mmathop{sum}limits_{jne i}{b}_{j}(s)frac{{a}_{i,j}}{{d}_{j}}{n}_{t}^{(j)}(s)+frac{1}{2}mu {sigma }_{mu }^{2}{{{Delta }}}_{s}left[{b}^{(i)}(s){n}_{t}^{(i)}(s)right]$$
(4)
Equation (4) is similar to Eq. (3), except that it incorporates an additional term corresponding to mutation processes and that the birth rate is trait-dependent. We show how Eq. (4) can be reduced to an equivalent two-habitat model under mean-field assumption. The mean-field approach differs slightly from the setting with no selection because vertices are labelled with Θi. Here we assume that vertices with similar habitat types have an equivalent position on the graph (see Supplementary Fig. 5 for a graphical representation), so that all vertices with habitat type I are characterised by the identical adaptive trait distribution that we denote by ({overline{n}}^{{{{{{{{bf{I}}}}}}}}}), and are associated with the birth rate ({b}^{{{{{{{{bf{I}}}}}}}}}(s)=b(1-p{(s-{theta }_{{{{{{{{bf{I}}}}}}}}})}^{2})). Let P(I, II) denote the proportion of edges connecting a vertex vi of type II to a vertex vj of type I, and let P(I) denote the proportion of vertices vi of type I. By further assuming that habitats are homogeneously distributed on the graph so that (P({{{{{{{bf{I}}}}}}}})=P({{{{{{{bf{II}}}}}}}})=frac{1}{2}), Eq. (4) transforms into$${partial }_{t}{overline{n}}_{t}^{{{{{{{{bf{I}}}}}}}}}(s)= ,{overline{n}}_{t}^{{{{{{{{bf{I}}}}}}}}}(s)left[{b}^{{{{{{{{bf{I}}}}}}}}}(s)(1-m)-frac{1}{K}{int}_{{{{{{{{mathcal{S}}}}}}}}}{overline{n}}_{t}^{{{{{{{{bf{I}}}}}}}}}({{{{{{{bf{s}}}}}}}})d{{{{{{{bf{s}}}}}}}}right]+frac{1}{2}mu {sigma }_{mu }^{2}({{{Delta }}}_{s}{b}^{{{{{{{{bf{I}}}}}}}}}{overline{n}}_{t}^{{{{{{{{bf{I}}}}}}}}})(s)\ +frac{m}{2},[(1-{r}_{{{Theta }}}){b}^{{{{{{{{bf{II}}}}}}}}}(s){overline{n}}_{t}^{{{{{{{{bf{II}}}}}}}}}(s)+(1+{r}_{{{Theta }}}){b}^{{{{{{{{bf{I}}}}}}}}}(s){overline{n}}_{t}^{{{{{{{{bf{I}}}}}}}}}(t)]$$
(5)
(see Methods), where we define$${r}_{{{Theta }}}=2left(P({{{{{{{bf{I}}}}}}}},{{{{{{{bf{I}}}}}}}})-P({{{{{{{bf{I}}}}}}}},{{{{{{{bf{II}}}}}}}})right)$$
(6)
as the habitat assortativity of the graph, which ranges from −1 to 1. When rΘ = − 1, all edges connect dissimilar habitat types (disassortative graph), while as rΘ tends towards 1 the graph is composed of two clusters of vertices with identical habitat types (assortative graph). Eq. (5) can be analysed with the theory of adaptive dynamics36,38,41, a mathematical framework that provides analytical insights by assuming a “trait substitution process”. Following this assumption, the mutation term in Eq. (5) is omitted and the phenotypic distribution results in a collection of discrete individual types that are gradually replaced by others until evolutionary stability is reached (see Methods and36,38,41 for details). By applying the theory of adaptive dynamics, we find a critical migration rate m⋆$${m}^{star }=frac{1}{(1-{r}_{{{Theta }}})}frac{4p{theta }^{2}}{(1+3p{theta }^{2})}$$
(7)
so that when m  > m⋆, a single type of individual exists with adaptive trait ({s}^{* }=left({theta }_{{{{{{{{bf{II}}}}}}}}}+{theta }_{{{{{{{{bf{I}}}}}}}}}right)/2=0) in the steady-state (see Methods for the derivation of Eq. (7)). In this case, adaptive differentiation QST,s is nil and the average population size is given by (overline{N}=bK{(1-ptheta )}^{2}). In contrast, when m = 0 and/or rΘ = 1, all individuals are locally well-adapted with trait Θi on vi, and it follows that the average population size is higher and equal to (overline{N}=bK), while adaptive differentiation is maximal and equal to ({Q}_{ST,s}={{{{{{{rm{Var}}}}}}}}({{Theta }})/left({{{{{{{rm{Var}}}}}}}}({{Theta }})+0right)=1). When 0  m⋆, implying that individuals become equally fit in all habitats. In this case, the isolation effect of heterogeneous selection is lost and QST,u reaches a similar level as in the setting with no selection for m  > m⋆ (Fig. 5a), although QST,u is slightly higher in the setting with heterogeneous selection due to lower population size ((overline{N}=bK(1-ptheta )) vs. (overline{N}=bK), see section above and Methods). This suggests that rΘ reinforces QST,u, as assortative graphs sustain higher levels of adaptive differentiation (Figs. 3 and 4). Simulations on the path graph with varying Θ-spatial distribution support this conclusion for high migration regimes, but show the opposite relationship under low migration regimes, where the habitat assortativity rΘ decreases QST,u (Fig. 5b). Assortative graphs are composed of large clusters of vertices with similar habitats, within which migrants can circulate without fitness losses. Local neutral trait distributions become more correlated within these clusters, resulting in a decline in QST,u for assortative graphs compared with disassortative graphs. Figure 5b therefore highlights the ambivalent effect of rΘ on QST,u. rΘ reinforces QST,u by favouring adaptive differentiation, but also decreases QST,u by decreasing population isolation within clusters of vertices with the same habitat type.We compare the effect of rΘ on QST,u to the effect of the topology metrics 〈l〉 and hd found in the setting with no selection using multivariate regression analysis on simulation results obtained for different graphs with varying Θ-spatial distribution (Fig. 5d for graphs with M = 7 vertices and Supplementary Fig. 7b for M = 9). The multivariate model explains the discrepancies in QST,u across the simulations for low and high migration regimes (see Supplementary Table 3 for details), and we find that rΘ, 〈l〉, and hd contribute similarly to neutral differentiation. Hence, the effects of rΘ and the topology metrics 〈l〉 and hd add up under heterogeneous selection. A change in sign of the standardized effect of rΘ on QST,s for low and high migration regimes verifies that the ambivalent effect of rΘ on QST,u found on the path graph holds for general graph ensembles. Simulations with trait-dependent competition and simulations on realistic graphs with a continuum of habitat types equally confirm the ambivalent effect of rΘ and further support the complementary effect of 〈l〉 and hd on QST,u (see Supplementary Fig. 8). 〈l〉 and hd therefore drive neutral differentiation with and without heterogeneous selection. rΘ becomes an additional determinant of neutral differentiation under heterogeneous selection. In contrast to the non-ambivalent, positive effect of habitat assortativity on adaptive differentiation, rΘ can amplify or depress neutral differentiation depending on the migration regime considered. More

The following section assesses our main results in terms of the growth in fenced areas over time relative to 1) types of protection, 2) administrative boundaries, and 3) other fences.Fencing relative to land governanceAcross the Greater Mara, a general growth in fenced areas can be observed throughout the 00 s but in particular over the last decade (Fig. 1). Based on satellite images, 35,067 ha were fenced in 1985, corresponding to c. 5%. In the following 25 years there was only an insignificant increase in fenced plots. However, from 2010, the number of fences suddenly grew rapidly, and in the following period (2015–2020) the fenced area increased even more radically, in an exponential manner (Fig. 2). For example, in 2015 there was 63,112 ha of fenced land; in 2016 this number rose to c. 75,176 ha, corresponding to a c. 20% annual increase. From 2010 to 2020, the ha fenced area increased by 170%. This corresponds to a roughly four times increase in the area enclosed by fences during the study period (1985–2020).Figure 2Conservative estimate of the fenced area of the entire Greater Mara, Kenya (1985–2020) expressed in hectares.Full size imageIn almost all regions, the number of fences continued to increase in 2019–20 (Fig. 2). The result is a total of 130,277 ha of fenced land in 2020, corresponding to 19% of the Greater Mara.Hence, there appears to be a building momentum in the expansion of fences in the Greater Mara: those regions that had many fences in 2016 ( > 1,000 ha) continue to experience an increase in the area enclosed by fences, with fences spreading almost everywhere in 2020 in particular. Those regions with the fewest fences in 2016 ( More

Ceballos, G. & Ehrlich, P. R. Mammal population losses and the extinction crisis. Science 296, 904–907 (2002).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Barnosky, A. D. et al. Has the Earth’s sixth mass extinction already arrived?. Nature 471, 51–57 (2011).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Ripple, W. J. et al. Status and ecological effects of the world’s largest carnivores. Science 343, 1241484 (2014).PubMed 
Article 
CAS 

Google Scholar 
Carbone, C. & Gittleman, J. L. A common rule for the scaling of carnivore density. Science 295, 2273–2276 (2014).ADS 
Article 

Google Scholar 
Durant, S. M. et al. The global decline of cheetah Acinonyx jubatus and what it means for conservation. Proc. Natl. Acad. Sci. 114, 528–533 (2017).CAS 
PubMed 
Article 

Google Scholar 
Jowkar, H. et al. Acinonyx jubatus ssp. venaticus. The IUCN Red List of Threatened Species 2008: e.T220A13035342. (2008).Khalatbari, L., Yusefi, G. H., Martínez-Freiría, F., Jowkar, H. & Brito, J. C. Availability of prey and natural habitats are related with temporal dynamics in range and habitat suitability for Asiatic Cheetah. Hystrix 29, 145–151 (2018).
Google Scholar 
Asadi, H. The Environmental Limitations and Future of the Asiatic Cheetah in Iran. (1997).CACP. Annual Report. (2014).Khalatbari, L., Jowkar, H., Yusefi, G. H., Brito, J. C. & Ostrowski, S. The current status of Asiatic cheetah in Iran. Cat News 66, 10–13 (2017).
Google Scholar 
Marker, L. L. et al. Ecology of free-ranging cheetahs. in Cheetahs: Biology and Conservation (eds. Marker, L. L., Boast, L. K. & Schmidt-Kuntzel, A.) 107–119 (Elsevier, 2017). doi:https://doi.org/10.1016/B978-0-12-804088-1.00008-3Hayward, M. W., Hofmeyr, M., O’Brian, J. & Kerley, G. I. H. Prey preferences of the cheetah (Acinonyx jubatus) (Felidae: Carnivora): morphological limitations or the need to capture rapidly consumable prey before kleptoparasites arrive?. J. Zool. 270, 615–627 (2006).Article 

Google Scholar 
Mills, M. G. L., Broomhall, L. S. & Toit, J. T. Cheetah Acinonyx jubatus feeding ecology in the Kruger National Park and a comparison across African savanna habitats: is the cheetah only a successful hunter on open grassland plains?. Wildlife Biol. 10, 177–186 (2004).Article 

Google Scholar 
Wachter, B., Jauernig, O. & Breitenmoser, U. Determination of prey hair in faeces of free-ranging Namibian cheetahs with a simple method. Cat News 44, 8–9 (2006).
Google Scholar 
Marker, L. L., Muntifering, J. R., Dickman, A. J., Mills, M. G. L. & Macdonald, D. W. Quantifying prey preferences of free-ranging Namibian cheetahs. South Afr. J. Wildl. Res. 33, 43–53 (2003).
Google Scholar 
Wacher, T. et al. Sahelo-Saharan Interest Group Wildlife Surveys, Part 4: Ahaggar Mountains, Algeria (March 2005). (2005).Thuo, D. et al. An insight into the prey spectra and livestock predation by cheetahs in Kenya using faecal DNA metabarcoding. Zoology 143, 125853 (2020).PubMed 
Article 

Google Scholar 
Broekhuis, F., Thuo, D. & Hayward, M. W. Feeding ecology of cheetahs in the Maasai Mara, Kenya and the potential for intra- and interspecific competition. J. Zool. 304, 65–72 (2018).Article 

Google Scholar 
Cooper, A. B., Pettorelli, N. & Durant, S. M. Large carnivore menus: factors affecting hunting decisions by cheetahs in the Serengeti. Anim. Behav. 73, 651–659 (2007).Article 

Google Scholar 
Mills, M. G. L. Living near the edge: A review of the ecological relationships between large carnivores in the arid Kalahari. African J. Wildl. Res. 45, 127–137 (2015).Article 

Google Scholar 
Rostro-García, S., Kamler, J. F. & Hunter, L. T. B. To kill, stay or flee: The effects of lions and landscape factors on habitat and kill site selection of cheetahs in South Africa. PLoS ONE 10, e0117743 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Laurenson, M. K. Behavioural costs and constraints of lactation in free-living cheetahs. Anim. Behav. 50, 815–826 (1995).Article 

Google Scholar 
Farhadinia, M. S. & Hemami, M.-R. Prey selection by the critically endangered Asiatic cheetah in central Iran. J. Nat. Hist. 44, 1239–1249 (2010).Article 

Google Scholar 
Farhadinia, M. S. et al. Feeding ecology of the Asiatic cheetah Acinonyx jubatus venaticus in low prey habitats in northeastern Iran: Implications for effective conservation. J. Arid Environ. 87, 206–211 (2012).ADS 
Article 

Google Scholar 
Zahedian, B. & Nezami, B. Cheetah (Acinonyx jubatus venaticus) (Felidae: Carnivora) feeding ecology in Central Plateau of Iran and effects of prey poor management. J. Wildl. Biodivers. 3, 22–30 (2019).
Google Scholar 
Zamani, N. et al. Predation of montane deserts ungulates by Asiatic cheetah Acinonyx jubatus venaticus in Central Iran. Folia Zool. 66, 50–57 (2017).Article 

Google Scholar 
Monterroso, P. et al. Factors affecting the (in)accuracy of mammalian mesocarnivore scat identification in South-western Europe. J. Zool. 289, 243–250 (2013).Article 

Google Scholar 
Morin, D. J. et al. Bias in carnivore diet analysis resulting from misclassification of predator scats based on field identification. Wildl. Soc. Bull. 40, 669–677 (2016).Article 

Google Scholar 
Caro, T. M. Cheetahs of the Serengeti Plains: Group Living in an Asocial Species (University of Chicago Press, 1994).
Google Scholar 
Floyd, T. J., Mech, L. D. & Jordan, P. A. Relating wolf scat content to prey consumed. J. Wildl. Manage. 42, 528–532 (1978).Article 

Google Scholar 
Jethva, B. D. & Jhala, Y. V. Computing biomass consumption from prey occurrences in Indian wolf scats. Zoo Biol. 23, 513–520 (2004).Article 

Google Scholar 
Pompanon, F. et al. Who is eating what: diet assessment using next generation sequencing. Mol. Ecol. 21, 1931–1950 (2012).CAS 
PubMed 
Article 

Google Scholar 
Taberlet, P., Coissac, E., Pompanon, F., Brochmann, C. & Willerslev, E. Towards next-generation biodiversity assessment using DNA metabarcoding. Mol. Ecol. 21, 2045–2050 (2012).CAS 
PubMed 
Article 

Google Scholar 
Mata, V. A. et al. How much is enough? Effects of technical and biological replication on metabarcoding dietary analysis. Mol. Ecol. 28, 165–175 (2019).CAS 
PubMed 
Article 

Google Scholar 
Shehzad, W. et al. Prey preference of Snow Leopard (Panthera uncia) in South Gobi Mongolia. PLoS ONE 7, e32104 (2012).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Monterroso, P. et al. Feeding ecological knowledge: the underutilised power of faecal DNA approaches for carnivore diet analysis. Mamm. Rev. 49, 97–112 (2019).Article 

Google Scholar 
Shehzad, W. et al. Carnivore diet analysis based on next-generation sequencing: Application to the leopard cat (Prionailurus bengalensis) in Pakistan. Mol. Ecol. 21, 1951–1965 (2012).CAS 
PubMed 
Article 

Google Scholar 
Thuo, D. et al. Food from faeces: Evaluating the efficacy of scat DNA metabarcoding in dietary analyses. PLoS ONE 14, e0225805 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Araujo, M. S., Bolnick, D. I. & Layman, C. A. The ecological causes of individual specialisation. Ecol. Lett. 14, 948–958 (2011).PubMed 
Article 

Google Scholar 
Balme, G. A., Roex, N., Rogan, M. S. & Hunter, L. T. B. Ecological opportunity drives individual dietary specialization in leopards. J. Anim. Ecol. 89, 589–600 (2020).PubMed 
Article 

Google Scholar 
Bolnick, D. I. et al. The ecology of individuals: incidence and implications of individual specialization. Am. Nat. 161, 1–28 (2003).MathSciNet 
PubMed 
Article 

Google Scholar 
Harrington, L. A., Harrington, A. L., Hughes, J., Stirling, D. & Macdonald, D. W. The accuracy of scat identification in distribution surveys: American mink, Neovison vison, in the northern highlands of Scotland. Eur. J. Wildl. Res. 56, 377–384 (2010).Article 

Google Scholar 
Weiskopf, S. R., Kachel, S. M. & McCarthy, K. P. What are snow leopards really eating? Identifying bias in food-habit studies. Wildl. Soc. Bull. 40, 233–240 (2016).Article 

Google Scholar 
Durant, S. M., Caro, T. M., Collins, D. A., Alawi, R. M. & Fitzgibbon, C. D. Migration patterns of Thomson’s gazelles and cheetahs on the Serengeti Plains. Afr. J. Ecol. 26, 257–268 (1988).Article 

Google Scholar 
Lindsey, P. A. et al. Minimum prey and area requirements of the vulnerable cheetah Acinonyx jubatus: implications for reintroduction and management of the species in South Africa. Oryx 45, 587–599 (2011).Article 

Google Scholar 
Farhadinia, M. S., Akbari, H., Eslami, M. & Adibi, M. A. A review of ecology and conservation status of Asiatic cheetah in Iran. Cat News Spec. Issue 18–26 (2016).Asadi, H. Some Observation on Hunting Behaviours of the Iranian Cheetah in Captivity. (1997).Heptner, V. G. & Sludskii, A. A. Mammals ofthe Soviet Union volume II part 2 Carnivora (hyaenas and cats). (Vysshaya Shkola Publishers, 1974).Ziaie, H. A Field Guide to the Mammals of Iran. (Iran Wildlife Center, 2008).Wilson, J. W. et al. Cheetahs, Acinonyx jubatus, balance turn capacity with pace when chasing prey. Biol. Lett. 9, 20130620 (2013).PubMed 
PubMed Central 
Article 

Google Scholar 
Grohé, C., Lee, B. & Flynn, J. J. Recent inner ear specialization for high-speed hunting in cheetahs. Sci. Rep. 8, 2301 (2018).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Cheraghi, F. et al. Inter-dependent movements of Asiatic Cheetahs Acinonyx jubatus venaticus and a Persian Leopard Panthera pardus saxicolor in a desert environment in Iran (Mammalia: Felidae). Zool. Middle East 65, 283–292 (2019).Article 

Google Scholar 
Ghoddousi, A., Soofi, M., Hamidi, A. K. & Lumetsberger, T. Assessing the role of livestock in big cat prey choice using spatiotemporal availability patterns. PLoS ONE 11, e0153439 (2016).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Khorozyan, I., Ghoddousi, A., Soofi, M. & Waltert, M. Big cats kill more livestock when wild prey reaches a minimum threshold. Biol. Conserv. 192, 268–275 (2015).Article 

Google Scholar 
Zeder, M. A. Domestication and early agriculture in the Mediterranean Basin: Origins, diffusion, and impact. Proc. Natl. Acad. Sci. 105, 11597–11604 (2008).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Daberger, M. Systematic prioritization of livestock grazing rights buyout in the last viable population of Asiatic cheetah (Acinonyx jubatus venaticus) in Iran. (Humboldt University Berlin, 2021).Wolf, C. & Ripple, W. J. Prey depletion as a threat to the world’s large carnivores. R. Soc. Open Sci. 3, 160252 (2016).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Melzheimer, J. et al. Communication hubs of an asocial cat are the source of a human—carnivore conflict and key to its solution. Proc. Natl. Acad. Sci. 117, 33325–33333 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Malakoutikhah, S., Fakheran, S., Tarkesh, M. & Senn, J. Assessing future distribution, suitability of corridors and efficiency of protected areas to conserve vulnerable ungulates under climate change. Divers. Distrib. 26, 1383–1396 (2020).Article 

Google Scholar 
Long, R. A., Donovan, T. M., Mackay, P., Zielinski, W. J. & Buzas, J. S. Comparing scat detection dogs, cameras, and hair snares for surveying carnivores. J. Wildl. Manage. 71, 2018–2025 (2007).Article 

Google Scholar 
Becker, M. S. et al. Using dogs to find cats: Detection dogs as a survey method for wide-ranging cheetah. J. Zool. 302, 184–192 (2017).Article 

Google Scholar 
Johnson, W. E. & O’Brien, S. J. Phylogenetic reconstruction of the Felidae using 16S rRNA and NADH-5 mitochondrial genes. J. Mol. Evol. 44, S98–S116 (1997).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Reese, E. M., Winters, M., Booth, R. K. & Wasser, S. K. Development of a mitochondrial DNA marker that distinguishes domestic dogs from Washington state gray wolves. Conserv. Genet. Resour. 12, 497–501 (2020).Article 

Google Scholar 
Ormerod, S. J. Applied issues with predators and predation: Editor’s introduction. J. Appl. Ecol. 39, 181–188 (2002).Article 

Google Scholar 
Boast, L. K., Good, K. & Klein, R. Translocation of problem predators: Is it an effective way to mitigate conflict between farmers and cheetahs Acinonyx jubatus in Botswana?. Oryx 50, 537–544 (2016).Article 

Google Scholar 
Darvish Sefat, A. A. Atlas of Protected Areas of Iran (University of Tehran, 2006).
Google Scholar 
Yusefi, G. H., Faizolahi, K., Darvish, J., Safi, K. & Brito, J. C. The species diversity, distribution, and conservation status of the terrestrial mammals of Iran. J. Mammal. 100, 55–71 (2019).Article 

Google Scholar 
Karami, M., Ghadirian, T. & Faizolahi, K. The Atlas of the Mammals of Iran. (Iran Department of the Environment, 2016).Abangah Consulting Engineer Company. Reconvene expanded Livestock Control Committee (LCC) in Touran and establish the LCC for Miandasht with participation of all stakeholders. (2017).Mills, M. G. L. & Hofer, H. Hyaenas. Status Survey and Conservation Action Plan. (IUCN/SSC Hyaena Specualist Group, 1998).Maudet, C., Luikart, G., Dubray, D., Von Hardenberg, A. & Taberlet, P. Low genotyping error rates in wild ungulate faeces sampled in winter. Mol. Ecol. Notes 4, 772–775 (2004).CAS 
Article 

Google Scholar 
Deagle, B. E., Kirkwood, R. & Jarman, S. N. Analysis of Australian fur seal diet by pyrosequencing prey DNA in faeces. Mol. Ecol. 18, 2022–2038 (2009).CAS 
PubMed 
Article 

Google Scholar 
Frantz, A. C. et al. Reliable microsatellite genotyping of the Eurasian badger (Meles meles) using faecal DNA. Mol. Ecol. 12, 1649–1661 (2003).CAS 
PubMed 
Article 

Google Scholar 
Boom, R. et al. Rapid and simple method for purification of nucleic acids. J. Clin. Microbiol. 28, 495–503 (1990).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Rosel, P. E. & Kocher, T. D. DNA-based identification of larval cod in stomach contents of predatory fishes. J. Exp. Mar. Bio. Ecol. 267, 75–88 (2002).Article 

Google Scholar 
Deagle, B. E. et al. Molecular scatology as a tool to study diet: analysis of prey DNA in scats from captive Steller sea lions. Mol. Ecol. 14, 1831–1842 (2005).CAS 
PubMed 
Article 

Google Scholar 
Riaz, T. et al. ecoPrimers: inference of new DNA barcode markers from whole genome sequence analysis. Nucleic Acids Res. 39, e145 (2011).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Luikart, G. et al. Multiple maternal origins and weak phylogeographic structure in domestic goats. Proc. Natl. Acad. Sci. 98, 5927–5932 (2001).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Menotti-Raymond, M. et al. A genetic linkage map of microsatellites in the domestic cat (Felis catus). Genomics 57, 9–23 (1999).CAS 
PubMed 
Article 

Google Scholar 
Charruau, P. et al. Phylogeography, genetic structure and population divergence time of cheetahs in Africa and Asia: Evidence for long-term geographic isolates. Mol. Ecol. 20, 706–724 (2011).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Driscoll, C. A., Menotti-Raymond, M., Nelson, G., Goldstein, D. & O’Brien, S. J. Genomic microsatellites as evolutionary chronometers: A test in wild cats. Genome Res. 12, 414–423 (2002).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Kotze, A., Ehlers, K., Cilliers, D. C. & Grobler, J. P. The power of resolution of microsatellite markers and assignment tests to determine the geographic origin of cheetah (Acinonyx jubatus) in Southern Africa. Mamm. Biol. 73, 457–462 (2008).Article 

Google Scholar 
Marker, L. L. et al. Molecular genetic insights on cheetah (Acinonyx jubatus) ecology and conservation in Namibia. J. Hered. 99, 2–13 (2008).CAS 
PubMed 
Article 

Google Scholar 
Taberlet, P. et al. Reliable genotyping of samples with very low DNA quantities using PCR. Nucleic Acids Res. 24, 3189–3194 (1996).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Egeter, B. et al. Challenges for assessing vertebrate diversity in turbid Saharan water-bodies using environmental DNA. Genome 61, 807–814 (2018).CAS 
PubMed 
Article 

Google Scholar 
Magoc, T. & Salzberg, S. L. FLASH: Fast length adjustment of short reads to improve genome assemblies. Bioinformatics 27, 2957–2963 (2011).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Rognes, T., Flouri, T., Nichols, B., Quince, C. & Mahé, F. VSEARCH: A versatile open source tool for metagenomics. PeerJ 4, e2584 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
Godinho, R. et al. Real-time assessment of hybridization between wolves and dogs: Combining noninvasive samples with ancestry informative markers. Mol. Ecol. Resour. 15, 317–328 (2015).CAS 
PubMed 
Article 

Google Scholar 
Valière, N. GIMLET: A computer program for analysing individual identification data. Mol. Ecol. 2, 377–379 (2002).
Google Scholar 
Wachter, B. et al. An advanced method to assess the diet of free-ranging large carnivores based on scats. PLoS ONE 7, e38066 (2012).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Breuer, T. Diet choice of large carnivores in northern Cameroon. Afr. J. Ecol. 43, 181–190 (2005).Article 

Google Scholar 
Wilson, M. F. J., O’Connell, B., Brown, C., Guinan, J. C. & Grehan, A. J. Multiscale terrain analysis of multibeam bathymetry data for habitat mapping on the continental slope. Mar. Geodesy 30, 2 (2007).Article 

Google Scholar 
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  More

Contents of Cd and Ca in Panax notogensing rootsThe Ca content of P. notoginseng roots increased significantly with the increase of lime application rates under the same concentration of oxalic acid sprayed on leaves (Table 2). Compared with no lime application, the Ca content was the highest increased by 212% under 3750 kg hm−2 lime without spraying oxalic acid. The content of Ca slightly increased with the increase of oxalic acid spraying concentrations under the same rate of lime application.Table 2 Effects of foliar spraying of oxalic acid on contents of Cd and Ca in roots of Panax notoginseng under Cd stress.Full size tableThe contents of Cd in roots ranged from 0.22 to 0.70 mg kg−1. The content of 2250 kg hm−2 Cd decreased greatly with the increase of lime application rates under the same spraying concentration of oxalic acid. Compared with the control, the root Cd contents decreased by 68.57% under the application of 2250 kg hm−2 lime and 0.1 mol L−1 oxalic acid spraying. The Cd contents of P. notoginseng roots decreased significantly with the increase of oxalic acid spraying concentrations under application of non-lime and 750 kg hm−2 lime. The root Cd contents decreased at first and then increased with the increase of oxalic acid concentrations under the application of 2250 kg hm−2 lime and 3750 kg hm−2 lime. In addition, the Bivariate analysis showed that the Ca content of P. notoginseng roots was significantly affected by lime (F = 82.84**), and the Cd content of P. notoginseng roots was significantly affected by lime (F = 74.99**) and oxalic acid (F = 7.72*).MDA contents and relative antioxidase activitiesThe content of MDA decreased greatly with the increase of the rates of lime application and oxalic acid spraying concentrations. There was no significant difference in the content of MDA in the roots of P. notoginseng with non-lime and 3750 kg hm−2 lime application. Under 750 kg hm−2, 2250 kg hm−2 lime application, the MDA content with 0.2 mol L−1 oxalic acid spraying concentration treatment decreased by 58.38% and 40.21% comparing with non-oxalic acid spraying application, respectively. The content of MDA (7.57 nmol g−1) was the lowest under 750 kg hm−2 lime application and 0.2 mol L−1 oxalic acid spraying treatment (Fig. 1).Figure 1Effects of foliar spraying of oxalic acid on contents of malondialdehyde in roots of Panax notoginseng under Cd stress. Notes The figure legend showed the spray concentration of oxalic acid (mol L−1), different lowercase letters indicate significant differences between treatments at the same lime application rate (P  Rb1  > R1. The contents of the three saponins had no significant difference with increase of the concentrations of oxalic acid spraying and no application of lime (Table 4).Table 4 Effects of foliar oxalate application on the percentages of three saponins in roots of Panax notoginseng under Cd stress.Full size tableThe contents of R1 with 0.2 mol L−1 oxalic acid spraying was significantly lower than that without oxalic acid spraying and rates of 750 or 3750 kg hm−2 lime application. Under the concentration of 0 or 0.1 mol L−1 oxalic acid spraying, there was no significant difference in contents of R1 with increase of rates of lime application. Under the concentration of 0.2 mol L−1 oxalic acid spraying, the contents of R1 with 3750 kg hm−2 lime was significantly lower 43.84% than that without lime application (Table 4).The contents of Rg1 increased at first and then decreased with the increase of oxalic acid spraying concentrations and 750 kg hm−2 lime application. Under the application rates of 2250 or 3750 kg hm−2 lime, the contents of Rg1 decreased with the increase of oxalic acid spraying concentration. With the same concentration of oxalic acid spraying, the Rg1 content increased at first and then decreased with the increase of lime application rates. Compared with the control, except that the Rg1 content with three concentrations of oxalic acid spraying and 750 kg hm−2 lime was higher than that of the control, the contents of Rg1 in the roots of P. notoginseng under other treatments was lower than that of the control. The Rg1 content was the highest with 750 kg hm−2 lime and 0.1 mol L−1 oxalic acid spraying treatment, which was higher 11.54% than that of the control (Table 4).The contents of Rb1 increased first and then decreased with the increase of oxalic acid spraying concentration and 2250 kg hm−2 lime application. The content of Rb1 with 0.1 mol L−1 oxalic acid spraying reached the maximum value of 3.46%, which was higher 74.75% than that without oxalic acid spraying treatment. Under other lime application treatments, there was no significant difference among different oxalic acid spraying concentrations. With 0.1 and 0.2 mol L−1 oxalic acid spraying treatments, the contents of Rb1 decreased at first and then decreased with the increase of lime application rates (Table 4).Contents of flavonoidsWith the same concentration of oxalic acid spraying, the content of flavonoids increased at first and then decreased with the increase of the amounts of lime application. There was no significant difference in the content of flavonoids under different concentrations of oxalic acid spraying without the application of lime or 3750 kg hm−2 lime. Under 750 and 2250 kg hm−2 lime application, the content of flavonoids increased at first and then decreased with the increase of the concentration of oxalic acid spraying. Under the treatment of 750 kg hm−2 application and 0.1 mol L−1 oxalic acid spraying, the content of flavonoids was the highest, which was 4.38 mg g−1, which was higher 18.38% than that of the same rate of lime application and without spraying oxalic acid. The content of flavonoids with 0.1 mol L−1 oxalic acid spraying treatment increased by 21.74% compared with that without oxalic acid spraying treatment and 2250 kg hm−2 lime application (Fig. 5).Figure 5Effects of foliar spraying of oxalate on the contents of flavonoids in roots of Panax notoginseng under Cd stress.Full size imageBivariate analysis showed that the content of soluble sugar in P. notoginseng root was significantly relationship with the amount of lime application and the concentration of oxalic acid spraying. The content of soluble protein in root was significantly relationship with lime application rates, both of lime and oxalic acid. The contents of free amino acid and proline in roots were significantly relationship with lime application rates, oxalic acid spraying concentrations, both of lime and oxalic acid (Table 5).Table 5 Variance analysis of the effects of oxalic acid, calcium and cadmium on the contents of multiple medicinal ingredients in the roots of Panax notoginseng (F value).Full size tableThe content of R1 in the root of P. notoginseng was significantly relationship with oxalic acid spraying concentrations, lime application rates, both of lime and oxalic acid. The content of flavonoids was significantly relationship with oxalic acid spraying concentrations, lime application rates. More

Pachauri, R. K. et al. Climate change 2014: synthesis report. Contribution of Working Groups I, II and III to the fifth assessment report of the Intergovernmental Panel on Climate Change. (IPCC, 2014).Thackeray, S. J. et al. Phenological sensitivity to climate across taxa and trophic levels. Nature 535, 241–245 (2016).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Walther, G.-R. et al. Ecological responses to recent climate change. Nature 416, 389–395 (2002).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Root, T. L., MacMynowski, D. P., Mastrandrea, M. D. & Schneider, S. H. Human-modified temperatures induce species changes: Joint attribution. Proc. Natl. Acad. Sci. USA 102, 7465–7469 (2005).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Menzel, A. et al. European phenological response to climate change matches the warming pattern. Glob. Change Biol. 12, 1969–1976 (2006).ADS 
Article 

Google Scholar 
Rosenzweig, C. et al. Attributing physical and biological impacts to anthropogenic climate change. Nature 453, 353–357 (2008).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Myneni, R. B., Keeling, C. D., Tucker, C. J., Asrar, G. & Nemani, R. R. Increased plant growth in the northern high latitudes from 1981 to 1991. Nature 386, 698–702 (1997).ADS 
CAS 
Article 

Google Scholar 
Walkovszky, A. Changes in phenology of the locust tree (Robinia pseudoacacia L.) in Hungary. Int. J. Biometeorol. 41, 155–160 (1998).ADS 
Article 

Google Scholar 
Crick, H. Q. P. & Sparks, T. H. Climate change related to egg-laying trends [8]. Nature 399, 423–424 (1999).ADS 
CAS 
Article 

Google Scholar 
Both, C., Van Asch, M., Bijlsma, R. G., Van Den Burg, A. B. & Visser, M. E. Climate change and unequal phenological changes across four trophic levels: Constraints or adaptations?. J. Anim. Ecol. 78, 73–83 (2009).PubMed 
Article 

Google Scholar 
Brown, J. L., Li, S. H. & Bhagabati, N. Long-term trend toward earlier breeding in an American bird: A response to global warming?. Proc. Natl. Acad. Sci. USA. 96, 5565–5569 (1999).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Miller-Rushing, A. J., Lloyd-Evans, T. L., Primack, R. B. & Satzinger, P. Bird migration times, climate change, and changing population sizes. Glob. Chang. Biol. 14, 1959–1972 (2008).ADS 
Article 

Google Scholar 
Hughes, L. Biological consequences of global warming: Is the signal already apparent?. Trends Ecol. Evol. 15, 56–61 (2000).CAS 
PubMed 
Article 

Google Scholar 
Kudo, G., Nishikawa, Y., Kasagi, T. & Kosuge, S. Does seed production of spring ephemerals decrease when spring comes early?. Ecol. Res. 19, 255–259 (2004).Article 

Google Scholar 
Doi, H., Gordo, O. & Katano, I. Heterogeneous intra-annual climatic changes drive different phenological responses at two trophic levels. Clim. Res. 36, 181–190 (2008).Article 

Google Scholar 
Both, C., Bouwhuis, S., Lessells, C. M. & Visser, M. E. Climate change and population declines in a long-distance migratory bird. Nature 441, 81–83 (2006).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Post, E. & Forchhammer, M. C. Climate change reduces reproductive success of an Arctic herbivore through trophic mismatch. Philos. Trans. R. Soc. B Biol. Sci. 363, 2369–2375 (2008).Article 

Google Scholar 
Thackeray, S. J. et al. Trophic level asynchrony in rates of phenological change for marine, freshwater and terrestrial environments. Glob. Chang. Biol. 16, 3304–3313 (2010).ADS 
Article 

Google Scholar 
Visser, M. E. & Both, C. Shifts in phenology due to global climate change: The need for a yardstick. Proc. R. Soc. B Biol. Sci. 272, 2561–2569 (2005).Article 

Google Scholar 
Hegland, S. J., Nielsen, A., Lázaro, A., Bjerknes, A. L. & Totland, Ø. How does climate warming affect plant-pollinator interactions?. Ecol. Lett. 12, 184–195 (2009).PubMed 
Article 

Google Scholar 
Brown, C. J. et al. Quantitative approaches in climate change ecology. Glob. Chang. Biol. 17, 3697–3713 (2011).ADS 
PubMed Central 
Article 

Google Scholar 
Parmesan, C., Duarte, C., Poloczanska, E., Richardson, A. J. & Singer, M. C. Overstretching attribution. Nat. Clim. Chang. 1, 2–4 (2011).ADS 
Article 

Google Scholar 
Damos, P. & Savopoulou-Soultani, M. Temperature-driven models for insect development and vital thermal requirements. Psyche (London) 2012, (2012).Osawa, T. et al. Climate-mediated population dynamics enhance distribution range expansion in a rice pest insect. Basic Appl. Ecol. 30, 41–51 (2018).Article 

Google Scholar 
Kiritani, K. Predicting impacts of global warming on population dynamics and distribution of arthropods in Japan. Popul. Ecol. 48, 5–12 (2006).Article 

Google Scholar 
Ozawa, A., Uchiyama, T. & Kasai, A. Estimating the day on which the numbers of adult tea spiny whiteflies (Aleurocanthus camelliae Kanmiya and Kasai) peaked and the number of generations produced in a year based on the effective cumulative temperature in tea fields. Annu. Rep. Kansai Plant Prot. Soc. 58, 57–64 (2016) (in Japanese).Article 

Google Scholar 
Ozawa, A., Saito, T. & Ikeda, F. Effect of host plant and temperature on reproduction of American Serpentine Leafminer, Liriomyza trifolii (Burgess). Japan. J. Appl. Entomol. Zool. 43, 41–48 (1999) (in Japanese).Article 

Google Scholar 
Baier, P., Pennerstorfer, J. & Schopf, A. PHENIPS—A comprehensive phenology model of Ips typographus (L.) (Col., Scolytinae) as a tool for hazard rating of bark beetle infestation. For. Ecol. Manage. 249, 171–186 (2007).Article 

Google Scholar 
Karuppaiah, V. & Sujayanad, G. K. Impact of climate change on population dynamics of insect pests. (2012).Lipper, L. et al. Climate-smart agriculture for food security. Nat. Clim. Chang. 4, 1068–1072 (2014).ADS 
Article 

Google Scholar 
Campbell, B. M. et al. Reducing risks to food security from climate change. Glob. Food Sec. 11, 34–43 (2016).Article 

Google Scholar 
Kiritani, K. The low development threshold temperature and the thermal constant in insects and mites in Japan (2nd edition). Bull. Nationai Instirute Agro-Environmental Sci. 31, 1–74 (2012) (in Japanese).
Google Scholar 
Nagasawa, A., Takahashi, A. & Higuchi, H. Host plant use for oviposition by Trigonotylus caelestialium (Hemiptera: Miridae) and Stenotus rubrovittatus (Hemiptera: Miridae). Appl. Entomol. Zool. 47, 331–339 (2012).Article 

Google Scholar 
Higuchi, H. Ecology and management of rice bugs causing pecky rice. Japan. J. Appl. Entomol. Zool. 54, 171–188 (2010) (in Japanese).Article 

Google Scholar 
Ohtomo, R. Occurrence and control of Stenotus rubrovittatus (Hemiptera: Miridae) in Touhoku area in Japan. Japan. J. Appl. Entomol. Zool. 57, 137–149 (2013) (in Japanese).Article 

Google Scholar 
Kiritani, K. The impact of global warming and land-use change on the pest status of rice and fruit bugs (Heteroptera) in Japan. Glob. Chang. Biol. 13, 1586–1595 (2007).ADS 
Article 

Google Scholar 
Nagasawa, A. & Higuchi, H. Suitability of poaceous plants for nymphal growth of the pecky rice bugs Trigonotylus caelestialium and Stenotus rubrovittatus (Hemiptera: Miridae) in Niigata, Japan. Appl. Entomol. Zool. 47, 421–427 (2012).Article 

Google Scholar 
Gordo, O. & Sanz, J. J. Temporal trends in phenology of the honey bee Apis mellifera (L.) and the small white Pieris rapae (L.) in the Iberian Peninsula (1952–2004). Ecol. Entomol. 31, 261–268 (2006).Article 

Google Scholar 
Sparks, T. H. & Yates, T. J. The effect of spring temperature on the appearance dates of British butterflies 1883–1993. Ecography (Cop.) 20, 368–374 (1997).Article 

Google Scholar 
Stefanescu, C., Penuelas, J. & Filella, I. Effects of climatic change on the phenology of butterflies in the northwest Mediterranean Basin. Glob. Chang. Biol. 9, 1494–1506 (2003).ADS 
Article 

Google Scholar 
Dell, D., Sparks, T. H. & Dennis, R. L. H. Climate change and the effect of increasing spring temperatures on emergence dates of the butterfly Apatura iris (Lepidoptera: Nymphalidae). Eur. J. Entomol. 102, 161–167 (2005).Article 

Google Scholar 
Tabuchi, K. et al. Rice bugs in the Tohoku Region: Their occurrence and damage from 2003 to 2013. Bull. Natl. Agric. Res. Cent. Tohoku Reg. 117, 63–115 (2015) (in Japanese).
Google Scholar 
Jolly, W. M. & Running, S. W. Effects of precipitation and soil water potential on drought deciduous phenology in the Kalahari. Glob. Chang. Biol. 10, 303–308 (2004).ADS 
Article 

Google Scholar 
Kriticos, D. J. et al. CliMond: Global high-resolution historical and future scenario climate surfaces for bioclimatic modelling. Methods Ecol. Evol. 3, 53–64 (2012).Article 

Google Scholar 
Takeda, A. & Shimizu, K. Characteristics of the damaged grain caused by the sorghum plant bug, Stenotus rubrovittatus (Hemiptera: Miridae) at different infection periods. Annu. Rep. Kanto-Tosan Plant Prot. Soc. 2009, 85–87 (2009) (in Japanese).
Google Scholar 
Ishimoto, M. Seasonal Prevalence of Occurrence of the Rice Leaf Bug, Trigonotylus caelestialium (Kirkaldy) (Heteroptera: Miridae) on Paddy Rice Plants. Jpn. J. Appl. Entomol. Zool. 48, 79–85 (2004) (in Japanese).
Katase, M., Shimizu, K., Siina, S., Hagiwara, K. & Iwai, H. Seasonal occurrence of rice bugs in the northern part of Chiba Prefecture. Annu. Rep. Kanto-Tosan Plant Prot. Soc. 2007, 99–104 (2007) (in Japanese).
Google Scholar 
Shintani, Y. Effect of seasonal variation in host–plant quality on the rice leaf bug, Trigonotylus caelestialium. Entomol. Exp. Appl. 133, 128–135 (2009).Article 

Google Scholar 
Takada, M. B., Yoshioka, A., Takagi, S., Iwabuchi, S. & Washitani, I. Multiple spatial scale factors affecting mirid bug abundance and damage level in organic rice paddies. Biol. Control 60, 169–174 (2012).Article 

Google Scholar 
Yoshioka, A., Takada, M. B. & Washitani, I. Landscape effects of a non-native grass facilitate source populations of a native generalist bug, Stenotus rubrovittatus, in a heterogeneous agricultural landscape. J. Insect Sci. 14, 1–14 (2014).Article 

Google Scholar 
Watanabe, T. & Higuchi, H. Recent occurrence and problem of rice bugs. Plant Prot. 60, 201–203 (2006) (in Japanese).
Google Scholar 
Seino, H. An estimation of distribution of meteorological elements using GIS and AMeDAS data. J. Agric. Meteorol. 48, 379–383 (1993) (in Japanese).Article 

Google Scholar 
Ishigooka, Y., Tsuneo, K., Nishimori, M., Hasegawa, T. & Ohno, H. Spatial characterization of recent hot summers in Japan with agro-climatic indices related to rice production. J. Agric. Meteorol. 67, 209–224 (2011).Article 

Google Scholar 
Yamasaki, K. et al. Intraspecific variations in life history traits of two pecky rice bug species from Japan: Mapping emergence dates and number of annual generations. Ecol. Evol. (2021).Sakagami, Y. & Korenaga, R. Triangle method—A simple method for the estimation of total effective temperature. Japan. J. Appl. Entomol. Zool. 25, 52–54 (1981) (in Japanese).Article 

Google Scholar 
Shigehisa, S. Seasonal changes in egg diapause induction and effects of photoperiod and temperature on egg diapause in the sorghum plant bug, Stenotus rubrovittatus (Matsumura) (Heteroptera: Miridae). Japan. J. Appl. Entomol. Zool. 52, 229–232 (2008) (in Japanese).Article 

Google Scholar 
Ishimoto, M. Oviposition of sorghum plant bug, Stenotus rubrovittatus (Matsumura) (Heteroptera: Miridae) on rice plants. Japan. J. Appl. Entomol. Zool. 55, 193–197 (2011) (in Japanese).Article 

Google Scholar 
Ogata, M. et al. Fecundity and longevity in adults of the sorghum plant bug, Stenotus rubrovittatus (Matsumura) (Heteroptera: Miridae) under laboratory conditions. Japan. J. Entomol. 13, 129–132 (2010) (in Japanese).
Google Scholar 
Niiyama, T. & Iitomi, A. Optimal control timing of the rice leaf bug, Trigonotylus caelestialium (Heteroptera: Miridae). Annu. Rep. Soc. Plant Prot. North Japan. (in Japanese) (2003).Takeda, A., Shimizu, K., Shiina, S., Hagiwara, K. & Katase, M. Seasonal prevalence of Stenotus rubrovittatus (Hemiptera: Miridae) in a gramineous weed field and a rice field. Annu. Rep. Kanto-Tosan Plant Prot. Soc. 2008, 97–102 (2008) (in Japanese).
Google Scholar 
Niiyama, T. Studies on the ecology of Trigonotylus caelestialium and establishment of pesticide control techniques. Annu. Rep. Akita Prefect. Agric. Inst. 49, 147–180 (2009) (in Japanese).
Google Scholar 
Sato, M. Latest infestation period of Miridae causing pecky rice damage in Aomori Prefecture. Annu. Rep. Soc. Plant Prot. North Japan 2014, 129–134 (2014) (in Japanese).
Google Scholar 
Johnson, J. B. & Omland, K. S. Model selection in ecology and evolution. Trends Ecol. Evol. 19, 101–108 (2004).PubMed 
Article 

Google Scholar 
Kuznetsova, A., Brockhoff, P. B. & Christensen, R. H. B. lmerTest package: tests in linear mixed effects models. J. Stat. Softw. 82, 1–26 (2017).Article 

Google Scholar  More

There are two phases involved in supporting the technical quality of the dataset: (i) the process used by Taronga when collecting data for each shark-bite incident, and (ii) the consistency modifications that we made during manuscript development.Phase oneFor each shark-bite incident, Taronga attempts to contact the most-relevant person involved in the event as possible — e.g., victim, victim’s family, witnesses. Those contacted are asked to complete a questionnaire with information relating to the shark bite. If applicable, questionnaires can also be completed by a fisheries officer in the relevant State and sent to Taronga. Taronga works closely with experts in each State’s fisheries department to validate information sourced in media reports. Each shark-bite case is unique, so the validation process varies depending on details specific to each incident. Forensic shark scientists within each State department are contacted after each shark bite to confirm details related to the incident. For example, validating species responsible for the bite often requires forensic analysis through expert examination of bite marks or artefacts21,22,23,24,25. When available, video footage is analysed for validation of information, such as confirmation of shark species and length.The database is cross-checked annually for the previous year with the International Shark-Attack File in Florida, USA, as well as with fisheries officers to ensure consistency. The database is cross-checked with the acknowledgement that there are discrepancies between versions due to differences in inclusion criteria. For example, the International Shark-Attack File includes bites where the victim was bitten aboard a boat, whereas the database we present here does not include bites aboard a boat.We acknowledge the limitations associated with this database, such as differences in reporting over time. For example, incidents might be reported more in recent decades due to technological advances making reporting more accessible or media publicising these events more widely18. There might also be reporting biases, for example, victims could be more likely to report a bite by a large, potentially dangerous shark (e.g., white, tiger, or bull shark) rather than a smaller, less-dangerous shark species (e.g., wobbegong shark). We also completed a quality assessment of the original database fields and redesigned the data acquisition and entry process (see Phase two) to allow exploration of shark-bite trends and patterns in Australia.Phase twoWe identified errors and inconsistencies in database fields. To avoid additional errors and inconsistencies, and to obtain a quality-controlled database, we redesigned the process for gaining and entering information into the database. This included creating a data descriptor (Supplementary File 2) used as a protocol to inform which questions to ask in the questionnaire. The data descriptor also directs the format of database entries by specifying information required in each field and by indicating the format of each entry (i.e., numeric, descriptive, or categorical). We manually inspected all previously entered data and adapted them to match the data descriptor. We checked each entry using the filter function in Microsoft Excel to identify any spelling and grammatical errors in the fields and ensure that all categories were grammatically identical. We standardised all metrics during this process (e.g., the data descriptor now stipulates that all length measurements should be recorded in metres).We validated and standardised the geographical locations of shark bites by converting all coordinates into decimal degrees using Microsoft Excel. We subsequently plotted all coordinates using the ggmap library27 in R (Version 4.0.2) (R Core Team 2020) (Fig. 2). We corrected any unusual coordinates (e.g., outside of Australian waters) and crosschecked them with site descriptions and states to ensure validity.Fig. 2Geographical locations of 1,196 shark bites in Australia. Each shark-bite incident is indicated by a red dot. (a) all shark-bite incidents; (b) bites most likely inflicted by bull sharks, (c) tiger sharks, and (d) white sharks. Two bite incidents that occurred at the Australian external territory, Cocos (Keeling) Islands, are not included on these maps. Background layers show elevation and major, perennial watercourses.Full size imageDuring the development of the data-descriptor protocol, we converted some previous descriptive columns into categorical columns. These columns included (but are not limited to) victim activity, attractant, injury location, injury severity, and weather condition. Converting these columns into categories facilitates analysis to investigate shark-bite patterns. For example, we converted victim activity into a categorical field to restrict answers to the following: snorkelling, motorised boating, unmotorised boating, boarding, swimming, standing, diving, fishing, or other, rather than allowing answers in any format. We used this information to create a time series to show the activity of shark-bite victims in Australia over time (Figs. 3 and 4). Shark bites have increased for boarders (including surfboarding, bodyboarding, kiteboarding, sailboarding, wakeboarding, and stand-up paddle boarding) over time, particularly since 1960 (Figs. 3b, 4). This is likely due to the increase in popularity of board sports, particularly surfing, since the 1960s28. This trend is likely not reflected in Fig. 3a because shark bites are unlikely to be classed as ‘provoked’ during board riding.Fig. 3Number of shark bites (black, dashed line) and proportion of activity done by the victim at the time of shark-bite incidents in Australia from 1900 to 2022. Panels represent shark bites in Australia that are; (a) provoked or (b) unprovoked. Boarding includes surfboarding, bodyboarding, kiteboarding, sailboarding, wakeboarding, and stand-up paddle boarding. Swimming includes snorkelling, spearfishing, freediving, body surfing, clinging to an object, falling into water, floating, or wading. Diving includes scuba-diving, hookah diving, or hard-hat diving. Fishing includes cleaning fish. No data for years 1908 and 1970.Full size imageFig. 4Number of shark bites (black, dashed line) and proportion of activity done by the victim at the time of shark-bite incidents in Australia from 1900 to 2022. Panels represent shark bites in Australia that are; (a) fatal or (b) non-fatal. Boarding includes surfboarding, bodyboarding, kiteboarding, sailboarding, wakeboarding, and stand-up paddle boarding. Swimming includes snorkelling, spearfishing, freediving, body surfing, clinging to an object, falling into water, floating, or wading. Diving includes scuba-diving, hookah diving, or hard-hat diving. Fishing includes cleaning fish. No data for years 1908 and 1970.Full size imageThe column representing a shark-bite victim’s recovery status also requires a categorical response, restricting answers to: fatal, injured, or uninjured (Fig. 5). Since 1900, the proportion of shark-bite-related fatalities has decreased (Fig. 5). This trend is also true for the three species most attributed to shark-bite-related fatalities, white (Carcharodon carcharias), tiger (Galeocerdo cuvier), and bull (Carcharhinus leucas) sharks. The decrease in shark-bite-related deaths is likely due to advancements in medical responses to shark-bite victims over time14 and better understanding among surfers about using tourniquets to stem bleeding following increased certification as first responders in workplace occupational health and safety requirements. Bites resulting in an uninjured victim includes interactions where the shark might have bitten the victim’s equipment (i.e., surfboard, bodyboard, kayak) rather than biting the person.Fig. 5Proportion of victim-recovery status (fatal = grey; red = injured; blue = uninjured) resulting from all unprovoked shark bites in Australia between 1900–2022. Blank years represent years without any reported occurrences. (a) all shark-bite incidents, (b) bites most likely inflicted by bull sharks, (c) tiger sharks, or (d) white sharks. No data for years 1908 and 1970.Full size imagePreviously, entries in the injury location column were descriptive. We converted the column to a categorical field restricting answers to the following: arm, hand, lower arm, upper arm, shoulder, neck, head, torso, leg, foot, calf, thigh, pelvic region, or other. During the analyses process, we further categorised these injury locations into four body areas (head, arm, torso, and leg) to assess how injury location affects recovery status (fatal or injured) (Fig. 6). Fatality most often occurred following shark bites to the torso (Fig. 6). This is likely due to the injuries to organs and major arteries resulting in blood loss, which is a leading cause of shark-bite fatalities29. This is the first time that the location of a shark bite on the body has been assessed relative to recovery status.Fig. 6Proportion of Australian shark bites resulting in either fatality or injury categorised by injury location on the victim’s body (left panel; 250 bites resulting in fatal injury, 723 bites resulting in non-fatal injury) and by species (right panel; bites by 201 tiger, 170 bull, 258 white, and 303 bites other sharks).Full size imageUnderstanding how the location of a shark-bite wound relates to victim recovery has value in informing the development of shark-bite mitigations. For example, the development of shark-bite-resistant wetsuits could potentially result in higher survival rates of the user if the fabric is concentrated around the torso region30,31. Redesigning data acquisition and entry process to allow for categorical columns permits these types of analyses.Some detail of a shark-bite incident might be lost by converting previously descriptive columns into categorical columns. We addressed this by complementing categorical columns with accompanying fields to retain details of the incident. For example, injury location and injury severity columns are both categorical and allow for user-friendly data analysis, whereas the injury description column is descriptive and provides added detail about the victim’s injuries if applicable. Individually, all three columns address certain aspects of the victim’s injuries, and together, all three columns comprehensively summarise injuries to the shark-bite victim.Our analysis of the Australian Shark-Incident Database suggests that tiger sharks are proportionally responsible for the most fatalities of all shark species in Australia (38% of all tiger shark bites result in fatality), followed by bull sharks (32% of all bull shark bites result in fatality), and white sharks (25% of all white shark bites result in fatality) (Fig. 6). We emphasise that these figures represent the overall percentage of bites resulting in fatality since 1791 and do not account for possible changes over time. These figures are also proportional to the number of bites by each respective species. In Australia, white sharks are responsible for the largest number of bites on humans (361 total) compared to tiger (229 total) and bull sharks (197 total). At the time of publication, white and tiger sharks were each responsible for 91 and 86 total fatalities on humans in Australia, respectively.There were 540 incidents in which time of day was recorded. We used these data to assess whether particular shark-bite incidents are more likely to occur at specific times of day (Fig. 7). To standardise reported 24-hour times, we took the location (latitude, longitude) information and reported time and date of each incident using the getSunlightTimes function in the suncalc library in R32 to calculate whether the incident occurred in one of four light-availability categories: dawn, day, dusk, or night. Shark bites occur mostly during the day, which likely reflects time of day when there are more ocean users present. However, there is a slightly higher proportion of bites at dusk for bull sharks compared to tiger and white sharks (Fig. 7). Identifying these trends can assist authorities in developing data-driven educational messaging as a shark-bite mitigation measure. This is important considering that enhanced education is the preferred mitigation measure scored by ocean users in New South Wales33.Fig. 7Period of day (corrected for local time) distribution of provoked and unprovoked shark bites in Australia by shark species from 1791 to 2022 (n = 540 bites from all species, 70 from bull, 59 from tiger, and 217 from white sharks).Full size imageThese examples demonstrate that the Australian Shark-Incident Database will be useful for scientists to analyse environmental, social, and biological related shark-bite patterns in Australia. Use of the newly developed data descriptor to standardise future applications and account for quality assurance and control will aid in keeping the database consistent for ease of analysis and interpretation. Ultimately, the publishing of this database will improve our understanding of shark-bite incidents in Australia and will equip us with the knowledge to aim to avoid or predict these events in the future. More