Philippa Kaur

Phylogenetic relationships inside the family Palaemonidae remain unresolved, despite being frequently discussed in recent publications9,10. Nevertheless, the last published study5 presented the main lineages of the family as well supported. Among those, the studied Pon-I group of predominantly free-living taxa is basal-positioned to the remaining genera of the former subfamily Pontoniinae, usually more specialised and associated with a wide range of hosts. The basal separation of the symbiotic genera led some authors to consider the assemblage, following Bruce22, to be a primitive group, or descendants of such7,23. Additionally, Gan et al.8 suggested that the taxa of the Pon-I group might be direct descendants of the ancestors of the former subfamily Pontoniinae, sharing the main plesiomorphies appearing frequently in former palaemonine taxa, e.g., the genera Brachycarpus, Leptocarpus, Macrobrachium, or Palaemon. The median process on the fourth thoracic sternite can be considered a plesiomorphic feature; indeed, it is a common symplesiomorphy of all Pon-I taxa, including Ischnopontonia and Anapontonia, for which the process was formerly reported as missing24 (its presence was confirmed in present examined specimens). In addition to that, the mandibular palp occurring in the genera Exoclimenella, Eupontonia, Palaemonella, and Vir25, or the presence of two arthrobranchs on the third maxilliped in Exoclimenella26, can also be considered plesiomorphic features.The Pon-I group’s internal relations have been unclear until now due to lower generic and species coverage in previous studies4,5,8. The present analysis based on a six-marker molecular dataset allows a deeper insight into the phylogenetic relationships of the study group involving all 11 currently recognised genera, and represented by 52 species, i.e. about 60% of the overall known species diversity of the group. The results provide a strong support for the monophyly and/or taxonomic validity of the current genera Exoclimenella, Anapontonia, Ischnopontonia, and suggest the monophyly of genera Harpilius and Philarius. Moreover, the results reveal non-monophyly of the most speciose genera Palaemonella and Cuapetes, as well as the species-poor Eupontonia. The genus Palaemonella was found to be paraphyletic owing to the nested species of the genera Eupontonia and Vir, which all share a common synapomorphy, the presence of the mandibular palp (mentioned above). Such conclusion was expressed also in the study of Chow et al.5.The present phylogenetic analysis confirmed that the genus Cuapetes is not monophyletic, as found to a lesser extent, in a few previous molecular studies 4,5,23. In this study, the genus Cuapetes was recovered in four separate genetic lineages. The type species C. nilandensis is nested in the Clade 1 along with C. johnsoni and C. seychellensis. This phylogenetic finding is in line with the study of Marin and Sinelnikov27, who indicated morphological differences between two of the above-mentioned species and most of the remaining species of the genus (respective of the present Clade 5, also covering C. grandis, the type species of the ex-genus Kemponia), and questioned the validity of the two latter generic names. The further genetic lineage is shown by the position of C. americanus nested in the eastern Pacific—Atlantic branch of the genus Palaemonella (Clade 3). This result is also supported by recent phylogenetic studies suggesting the different systematic positions of this species4,5,10. Due to the lack of the mandibular palp, the species had been properly, but evidently incorrectly, assigned to the genus Cuapetes. The fourth genetic lineage is shown by a separate position of C. darwiniensis in the Clade 4 as the sister species of Madangella altirostris.The remaining majority of the Cuapetes species (Clade 5) are heterogeneous due to comprising also representatives of the genus Periclimenella. Ďuriš and Bruce26 hypothesised, based on morphological traits (mainly the unique shape of the first pereiopod chelae and the distinctly asymmetrical and specific second pereiopods), that the genera Exoclimenella and Periclimenella are closely related. Nevertheless, the present study revealed Periclimenella as a part of the genus Cuapetes. This result was previously supported in the molecular study by Horká et al.4 and weakly supported by Kou et al.23.Fossil records of palaemonid shrimps are rare due to their aquatic habit and poorly calcified exoskeletons. Only a few palaemonid representatives are known compared to many extant taxa; the oldest fossil records contain only genera from the previous subfamily Palaemoninae from the Lower Cretaceous (middle Albian, 100 Myr)28. For this reason, we used the known mutation rate of mitochondrial gene (16S rRNA) for dating rather than fossil records.The present inferred phylogeny and ancestral analysis indicate multiple formations of primary symbioses within the clades dominated by free-living relatives, as shown by previous molecular analyses4,5. Our results revealed eight independent lineages within the Pon-I group that evolved from free-living ancestors (Fig. 3). Free-living palaemonids (Exoclimenella, Palaemonella, Cuapetes; Fig. 2) are characterised by an elongate body shape with a dentate rostrum, slender, long, a/symmetrical chelipeds and slender ambulatory pereiopods with simple dactyli. Their carapace might bear the full complement of teeth (i.e., supraorbital, antennal, hepatic, epigastric)25. Primary symbiotic forms do not fundamentally differ morphologically from free-living ancestors. Their adaptations to the host affiliation have mainly manifested by changes in body shape, colouration, and the reduction of carapace ornamentation. Their hosts belong to different invertebrate phyla, including Cnidaria (mainly Scleractinia and Antipatharia22) and Echinodermata (Crinoidea29) in ectosymbiotic forms, but also to spoon worms (Echiura), burrowing Crustacea (alpheid shrimps), and/or gobiid fishes15, in inquilinistic forms.While scleractinian corals were hypothesised as the primary hosts of palaemonid shrimp commensalism7, our results revealed the antipatharian association as possibly the earlier one among the Pon-I shrimps. That association was established via a single speciation act at approximately 43 Myr (Eocene), specifically with the ancestor of the recent Cuapetes nilandensis (Clade 1). Except a small body size, this species does not show specific morphological adaptations to antipatharian association. The possibly oldest lineage associated with the scleractinian corals forms a common multigeneric composition of Anapontonia, Ischnopontonia, Harpilius and Philarius (Clade 4), which was established at approximately 38.2 Myr (Eocene). The genera share some homoplasic adaptations with ectosymbioses, such as strongly hooked dactyli of the ambulatory pereiopods adapted to climbing on coral colonies. An extremely compressed body and similar tail fan structure of the genera Ischnopontonia (Fig. 1H) and Anapontonia (Fig. 1D) are adaptations to life in narrow spaces amongst corallites of the oculinid coral Galaxea24,30; the intercorallite channels might be temporarily fully covered by tentacles of exposed polyps. This lifestyle was thus termed ‘semi-endosymbiosis’ by Horká et al.4, as potential evolutionary precursors of the true endosymbioses. In contrast, the genera Philarius and Harpilius have depressed bodies and associate exclusively as regular ectosymbionts with scleractinian corals, mainly of the genera Acropora and Pocillopora22.A further multispecies symbiotic lineage is represented by the genus Vir (Clade 3), whose origin is dated to approximately 21.1 Myr (Miocene). All species of this genus live in associations mainly with the acroporid, pocilloporid and euphylliid genera of scleractinian corals31,32. The adaptation to their symbiotic lifestyle is expressed in the loss of the hepatic tooth, partial or full reduction of ambulatory propodal spines, and cryptic colouration, including transparency of the body and appendages31,33 (Fig. 1J). Subsequent scleractinian-associated lineages are represented by separate species that appeared in the Miocene (21.9–10.1 Myr), namely: Eupontonia oahu, Cuapetes amymone, and C. kororensis, which live in association with Pocillopora, Acropora, and Heliofungia, and show only minor adaptations to their symbiotic habits, e.g. loss of the hepatic tooth, dense distal setae on the walking propodi, or extremely slender chelae and a specific cryptic colouration, respectively22,34,35.A single crinoid-associated species, Palaemonella pottsi (Clade 3), represents the only case of the switch from a free-living lifestyle to the association with echinoderms in the present study group; it originated at approximately 10.4 Myr (Miocene). Retaining the body shape typical for Palaemonella12, the species also does not show any noticeable morphological adaptation to such a host; its affiliation with the symbiotic life is, however, clearly observed in the deep-red to black cryptic colouration36.In Palaemonella aliska (Fig. 1E) and Eupontonia nudirostris (Clade 3), a pair of sister-positioned species in the present analyses (Figs. 2, 3), the ability to co-habit with burrowing animals (e.g., alpheids, gobiid fish, or echiurids) had developed. Their type of symbiosis, inquilinism, formed at approximately 14.8 Myr (Miocene). The reduction of the rostrum length, depressed body, stout main chelae in both, and full lack of the epigastric and hepatic teeth in the latter species15,25, were evidently due to that mode of life. Inquilinism is best known in the family Alpheidae, in which multiple genera associate with a variety of burrowing animals37. In the family Palaemonidae, inquilinism developed only in the Pon-I group, including Palaemonella shirakawai (not analysed here)14.As evident from the present and previously published reports4,5,7,8,10, the life history of the Pon-I group was largely shaped by coevolution with coral reefs. The coral reefs were deeply impacted by the K–T mass extinction at the end of the Cretaceous, which was one of the most destructive events in the Phanerozoic38. However, coral reefs recovered and became increasingly abundant in the Eocene39. This also matches the time of either the origin of host associations, or a wider species radiation of the Pon-I group. The first fossil records of the main coral hosts of the present shrimps are dated after the K-T extinction during the Paleogene (e.g., Euphyllia 66.0–61.6 Myr, Acropora 59.2–56.0 Myr, Galaxea and Pocillopora 56–33.9 Myr40).The biogeographic history suggested by S-DIVA analysis points to some dispersal and vicariant events shaping the current pattern of the Pon-I group’s distribution. This reconstruction (Fig. 4) estimates the present-day IWP region within the former Paleo-Tethys Ocean as the most likely ancestral area of the present study group, which originated ~ 91.6 Myr (Late to Early Cretaceous). The present shrimp group had radiated across the entire IWP region and subsequently expanded into the Atlantic Ocean. We assume that the spread of the group took place in the following sequence of events: (1) dispersal of Palaemonella spp. from the IWP into the eastern Pacific in the Paleocene (∼ 55.2 Myr; P. asymmetrica and P. holmesi); (2) dispersal into the western Atlantic (2 spp., complex of “Cuapetes” americanus) via the eastern Pacific and vicariance event separating the IWP at Eocene (∼ 46.2 Myr). It was the time after the formation of the Eastern Pacific Barrier (EPB), which was considered the largest extension of the open ocean (ca. 5000 km), that separated the IWP area from the eastern Pacific17; (3) the another vicariance event, separating the western Atlantic populations from those of the eastern Pacific in the Oligocene (∼ 30.9 Myr), i.e., before the closure of the Isthmus of Panama, followed by a dispersion of P. atlantica into the eastern Atlantic in the Miocene (∼ 21.6 Myr). The exact time of the formation of the Isthmus of Panama, which separated the Atlantic from the eastern Pacific and remained isolated from the central Pacific by the EPB, still remains questionable. Bacon et al.18 assume that the initial land bridge formed at approximately 23 Myr, and the final closure of the Isthmus of Panama formed between 10 and 6 Myr. Montes et al.19 presupposed the earlier formation of the barrier at ∼ 14 Myr, whereas O’Dea et al.20 concluded that the potential gene flow continued between the Pacific and Atlantic subpopulations of marine organisms until at least ∼ 2.8 Myr.The eastern Pacific Cuapetes canariensis closely related to IWP Cuapetes spp., has been recently described by Fransen et al.41, from the Canary Islands. This could indicate alternative dispersal pathways into the Atlantic, as suggested by recent studies17,42. The Tethys seaway allowed natural dispersion between the Atlantic and Indian Oceans across the region of the Mediterranean Sea. The closure of this interoceanic seaway at approximately 14 Myr (18–12 Myr) was caused by intense tectonic activity in the Near East17. Since the closure of that seaway, remaining possible dispersal to the Atlantic has been limited to the warm-water corridor around the southern tip of Africa, however curtailed by the cold Benguela Current upwelling from the Late Pliocene43. More

Study area and soil propertiesA field experiment was established in May 2012 at Shenyang Agro-Ecological Station (41°31′N, 123°22′E) of the Institute of Applied Ecology, Chinese Academy of Sciences, Northeast China. This region has a warm-temperate continental monsoon climate. The mean annual air temperature and annual precipitation are 7.5 °C and 680 mm, respectively. The soil is classified as Luvisol (FAO classification). The soil properties of the topsoil layer (0–20 cm) at the start of the experiment are as follows: SOC = 9.0 g kg−1, available NH4+–N = 1.18 mg kg−1; available NO3−–N = 9.04 mg kg−1; Olsen-P = 38.50 mg kg−1, available K = 97.90 mg kg−1, bulk density = 1.25 g cm−3, and pH = 5.8. The determination method of soil was shown in “Soil analysis” section.Field experimentThree treatments were established in this experiment: (1) mineral fertilizers (NPK); (2) pig manure incorporation at a local conventional AM application rate of 15 Mg ha−1 yr−1 (NPKM, 126 kg N ha−1 on dry weight); and (3) NPKM plus DMPP (3,4-Dimethylpyrazole phosphate) incorporation at a rate of 0.5% of applied urea (2.39 kg ha−1, 220 kg N/the N content of urea (0.46) × 0.5%) (NPKI + M). The treatments were applied following a randomized design across three replicate field plots (4 m × 5 m). Plots of different treatments remained unchanged in the same locations for 4 years. Each year, the composted pig manure (213 g C kg−1 and 22 g N kg−1 based on dry weight on average, characteristics of pig manure was listed in Table S1) was broadcasted evenly onto the plots a few days before maize planting, and ploughed to a depth of 20 cm by machine (TG4, Huaxing, China). For the respective treatments, urea (220 kg N ha−1 yr−1), calcium superphosphate (110 kg P2O5 ha−1 yr−1), and potassium chloride (110 kg K2O ha−1 yr−1) were applied on the same day as maize (Zea mays L.) was planted. The urea and inhibitor were fully mixed before application.Maize (cultivar was Fuyou #9) was planted on 3rd May 2012, 3rd May 2013, 6th May 2014, and 10th May 2015, at a spacing of 37 cm and 60 cm between rows. No irrigation was applied throughout the experimental period. Maize was harvested on 13th September 2012, 29th September 2013, 29th September 2014, and 29th September 2015, respectively. At harvest, maize yield and aboveground biomass yield were measured by harvesting all plants (20 m2) in each plot. The straw and grain were removed after each harvest and the soil with about 5 cm maize stem was ploughed to a depth of approximately 20 cm in April each year.Each cropping cycle, therefore, consisted of periods of maize (from May to September) and fallow (from October to April) of the following year.The precipitation and air temperature data were acquired from the meteorological station of the Shenyang Agro-Ecological Station. The precipitation during the 2012/2013, 2013/2014, 2014/2015, and 2015/2016 periods were 911.9 mm, 621.7 mm, 485.7 mm, and 585.3 mm, respectively (Fig. 1). 72.3%, 75.5%, 66.5%, and 73.0% of these annual precipitations occurred during maize-growing period, respectively. The mean annual air temperatures in these years were 7.7 °C (− 21.2 to 27.5 °C), 8.1 °C (− 22.7 to 28.3 °C), 9.5 °C (− 21.7 to 28.2 °C) and 9.3 °C (− 17.1 to 27.0 °C), respectively. The soil temperature at a depth of 5 cm varied between − 14 and 35 °C during the four-year period (Fig. 2b). The change trend of soil surface temperature was the same as that of soil temperature at 5 cm depth (Fig. 2a). The mean soil WFPS (0–15 cm) varied between 15 and 73% (Fig. 2c).Figure 1Precipitation and daily mean air temperature during four annual cycles from May 2012 to April 2016 in the experimental field.Full size imageFigure 2Seasonal variations in soil temperature (at soil surface and 5 cm soil depth) and WFPS% at 0–15 cm depth from May 2012 to April 2016.Full size imageGas sampling and analysisThe gas was sampled between 3rd May 2012 and 14th April 2016 using a static closed chamber system as described by Dong et al.16. Briefly, a stainless-steel chamber base (56 cm length × 28 cm width) was inserted into the soil of each plot to a depth of approximately 10 cm, with its long edge perpendicular to the rows of maize. The top chamber (56 cm length × 28 cm width × 20 cm height) was also made of stainless steel. Gas samples were obtained using a syringe 0, 20, and 40 min after the chambers had been closed between 9:00 am and 11:00 am on each sampling day. Gas samples were collected every 2‒6 days and every 7‒15 days during the growing seasons and non-growing seasons, respectively. The first gas sampling time was on day 1, day 3, day 1, and day 3 after maize planting each year. The N2O concentrations in gas samples were quantified using a gas chromatograph (Agilent 7890A, Shanghai, China) with an electron capture detector.Soil analysisThe soil temperature and volumetric water content (SVWC) were measured at depth of 0–15 cm using a bent stem thermometer and a time-domain reflectometry (Zhongtian Devices Co. Ltd, China), respectively. SVWC was converted to soil water-filled pore space (WFPS) using the following equation:$${text{WFPS}} = {text{SVWC}}/(1{-}{text{BD}}/{text{particle}},{text{density}}),$$
(1)
where BD is soil bulk density (g cm−3). Particle density was assumed to be 2.65 g cm−3.Soil samples from the 0–20 cm layer were collected in each plot in April 2012 (before sowing) and October 2015 (maize harvest) using a 5 cm diameter stainless steel soil sampler. The five soil samples collected from different locations in each plot were mixed thoroughly. Visible roots were removed by hand and the samples were air-dried and sieved using a 0.15 mm sieve. SOC was then quantified using an elemental analyzer (Vario EL III, Elementar, Germany). Soil available NH4+–N and NO3−–N were extracted with 2 M KCl and measured colorimetrically using a continuous flow injection analyzer (Futura, Alliance, France)17. Soil Olsen-P was extracted with NaHCO3 and colorimetrically measured using a spectrophotometer (Lambda 2, PerkinElmer, USA). Soil available K was extracted by 1 M CH3COONH4 and analyzed with a flame photometer (FP640, Jingmi, China). Soil pH was determined with deionized water (1:2.5) and analyzed using a pH meter (PHS-3C, LeiCi, China) with a glass electrode.DNA extraction and real-time quantitative PCRThe soil samples for measuring the abundance of nitrification and denitrification functional genes were collected on May 20, 2015. Soil DNA was extracted with the soil DNA extracted kits (EZNA soil DNA Kit; Omega Bio-Tek Inc., U.S.A.). The copy numbers of nitrification and denitrification functional genes were determined by q-PCR with the Roche LightCyler® 96 (Roche, Switzerland). Additional details about the primers and amplification procedure can be found in Dong et al.16.Data analysisThe N2O flux (μg N2O–N m−2 h−1) is calculated based on the increase of N2O concentration per unit chamber area for a specific time interval18 as follows:$${text{F}} = 273/left( {273 + {text{T}}} right) times {text{M}}/22.4 times {text{H}} times {text{dc}}/{text{dt}} times 1000$$
(2)
where F (μg N2O–N m−2 h−1) is the N2O flux, T (◦C) is the air temperature in the chamber, M (g N2O–N mol−1) is the molecular weight of N2O–N, 22.4 (L mol−1) is the molecular volume of the gas at 101.325 kPa and 273 K, H (m) is the chamber height, dc/dt (ppb h−1) is the rate of change in the N2O concentration in the chamber.Cumulative N2O emissions were calculated as follows:$${text{Cumulative}},{text{emission}} = mathop sum limits_{{{text{i}} = 1}}^{{text{n}}} frac{{({text{F}}_{{text{i}}} + {text{F}}_{i + 1} )}}{2} times ({text{t}}_{{{text{i}} + 1}} – {text{t}}_{{text{i}}} ) times 24$$
(3)
where F is the N2O emission flux (μg N2O–N m−2 h−1), i is the ith measurement, (ti+1 − ti) is the number of days between two adjacent measurements, and n is the total number of the measurements. Annual N2O emissions were calculated between the fertilization dates of each successive year.The SOC stock (Mg ha−1) in the topsoil was calculated as:$${text{C}}_{{{text{stock}}}} = {text{SOC}} times {text{BD}} times {text{D}} times 10,$$
(4)
where BD is soil bulk density (g cm−3), D is the depth of the topsoil (0.2 m).The topsoil SOC sequestration rate (SOCSR) (Mg ha−1 yr−1) was estimated using the following equation:$${text{SOCSR}} = left( {{text{C}}_{{{text{stock2015}}}} – {text{C}}_{{{text{stock2012}}}} } right) times {text{t}}^{ – 1} ,$$
(5)
where Cstock2015 and Cstock2012 are the SOC stocks in 2015 and 2012, respectively, and t is the duration of the experiment (years).Statistical analyses were performed using SPSS 13.0 (SPSS, Chicago, USA). The differences in cumulative N2O emissions and maize yields within a year, and other factors among treatments were assessed using one-way Analysis of Variance (ANOVA) with least significant difference post-hoc tests and a 95% confidence limit. The effects of different treatments, years, and their interactions on N2O emission, maize yield and aboveground biomass were examined using one-way repeated measures ANOVA. Pearson correlation analysis was used to analyze the relationships between cumulative N2O emissions and precipitation (N = 12 (three data each year, four years)), as well as N2O flux and soil available nitrogen content.
Statements of research involving plantsIt is stated that the current research on the plants comply with the relevant institutional, national, and international guidelines and legislation. It is also stated that the appropriate permissions have been taken wherever necessary, for collection of plant or seed specimens. It is also stated that the authors comply with the ‘IUCN Policy Statement on Research Involving Species at Risk of Extinction’ and the ‘Convention on the Trade in Endangered Species of Wild Fauna and Flora’. More

Data sources and pre-processingEach of the predictor variables used in our analysis (Table 1), as well as the dependent variable (fire hotspots) underwent pre-processing to transform the data into a format suitable to be passed to our CNN model for prediction. Here we briefly outline these processes and describe the method of generating a training and validation data set for model development. For further details about each predictor variable pre-processing, see Horton et al. (2021).Table 1 Model input data sources, citation, original resolution, and date ranges.Full size tableFire hotspotsWe used both Moderate Resolution Imaging Spectroradiometer (MODIS) and Visible Infrared Imaging Radiometer Suite (VIIRS) fire hotspot data as the dependent variable for use in our model development. As fire hotspots do not give precise locations, but rather indicate that a fire hotspot occurred within a grid cell of the size of the dataset (MODIS 1 km, VIIRS 375 m), we represented each fire hotspot as a 500 m buffered area around the centre point of each grid square identified. We used all fire hotspot occurrences with a confidence rating >50%.LandcoverWe use a collection of historic land cover maps generated by the Ministry of Forestry Indonesia from 1996 to 2016 at 2–3 year intervals38. Before use, we re-designated the land cover map classifications to reduce the number from 25 to just 8 (supplementary Table S2), which are ‘Primary and secondary dry forest’, ‘Swamp forest, ‘Swamp scrubland’, ‘Scrubland, Transition, and bare land’, ‘Riceland’, ‘Plantation’, ‘Settlements’, ‘water, and Cloud’.In addition to these 8 land cover classifications, we also derived a forest clearance index, which identifies areas cleared of forest and assigns an index value that is large negative (−10) immediately after clearing and degrades back towards 0 as time since clearing increases yearly. Areas that are re-forested are assigned large positive values (10) that degrade towards 0 yearly as time since afforestation increase25.Vegetation indicesAll vegetation indices were taken as pre-fire season 3-month averages from May to July. In addition to the original MODIS ET, PET, NDVI, and EVI products, we also included ‘normalised’ variables, whereby each vegetation index was expressed as the ratio of the same index taken at a reference site. The reference site was an area of dense primary forest outside of the EMRP area.Proximity to anthropogenic factorsThe distance to roads and settlement rasters were derived from OpenStreetMap data as the Euclidean distance to nearest feature in 250 m resolution. The same was done for all water bodies, which were then classified by hand into either canals or rivers. These features are taken as those shown in 2015 for all years, and therefore may misrepresent earlier years. However, the majority of canal development in the region took place between 1996 and 1998 and so should not differ dramatically from this date onwards.Oceanic Niño Index (ONI)We use a single value for the entire study area taken as the three-month average for the early fire season each year (July–September).Number of cloud daysUsing the state_1km band in the daily MODIS terra product (MOD09GA version 6), which classifies each pixel as either ‘no cloud’, ‘cloud’, ‘mixed’, or ‘unknown’, we counted the number of ‘cloud’ or ‘mixed’ designations for each pixel for the pre-fire season period May–July.Cross year normalisationAll predictor variables are normalised to be represented between 0 and 1 as the range between the minimum and maximum values for each variable that occur across all years, such that:$${V}_{{{{{{rm{norm}}}}}}}=frac{V-{V}_{{min }}}{{V}_{{max }}-{V}_{{min }}}$$where ({V}_{{{{{{rm{norm}}}}}}}) is the normalised version of the predictor variable (V), ({V}_{{max }}) is the maximum value within the training dataset across all years (2002–2019), and ({V}_{{min }}) is the minimum value within the training dataset across all years.Training and validation dataset assemblyOnce pre-processed, all predictor variable rasters were resampled to the same dimensions (with a resolution of 0.002 degrees in the WGS84 co-ordinate system) and stacked yearly, so that each year (2002–2019) comprised of a 31 feature maps input as a raster stack, with each feature map representing a different predictor variable. Each yearly stack was then split into tiles matching the input dimensions of the CNN model. Our final model was built to take an input size of 32 × 32 pixels (raster cells). Therefore, each yearly raster stack was split into many 32 × 32 × 31 raster stack tiles that span the defined study area. These were then converted to 3D arrays holding the values of all predictor variables for each raster stack tile.The same process was repeated for the yearly fire hotspot rasters used as the dependent variable in building our model. Each year was split into 32 × 32 × 1 tiles across the study area, and then converted to 3D arrays, each of which pairs with one predictor variable array.The 3D predictor variable arrays (dimensions: 32 × 32 × 31) were then stacked into one large 4D array containing all these individual tiles across all years (dimensions: W × 32 × 32 × 31, where W is a large value). The same was done with the 3D dependent variable arrays (dimension: 32 × 32 × 1), preserving the order so that each element in this large 4D array (dimensions: W × 32 × 32 × 1) matches with its counterpart in the predictor variable array.The order of this large 4D training data array was then randomised along the first dimension to avoid bias in passing to the CNN training algorithm, but the randomised re-ordering was repeated with the dependent variable array so as to preserve the elementwise pairing for cross-validation.Model development and applicationFire prediction requires the combination of spatial and temporal indicators to generate a probabilistic output for each location within a given study area. There is a need to preserve a certain level of proximity information, as the location of variables in relation to one another may have a substantial impact on the results. For example, a patch of secondary forest that is immediately adjacent to an area recently deforested may have a significantly higher probability of fire occurrence than an area surrounded entirely by primary forest.CNNs retain spatial features by employing a moving window of reference, known as a kernel, over the input image that captures these proximity relationships within the model structure. For this reason, CNNs are often used for image classification problems, and is an ideal model configuration for the problem of fire prediction across an area. Therefore, we have developed a CNN binary classification model using the Keras API package39 that builds on the TensorFlow machine learning platform40.Model structureCNN models typically apply a combination of kernel layers and dense layers that perform a series of transformations on the multi-channel input to either reduce it down to a single value, or to output an image the same width and height as the input with a single channel. These classification models can either assign a single value (binary classifier), or return one of many possible classifications.Kernels act on a subsection of the input stack (31 feature maps), assigning weights according to each cell’s position within the subsection to transform and combine the values into a new format to pass forward. As the kernel is applied to all subsections of the input stack, it transforms them to the new format, and builds a reconstituted image with dimensions that usually differ from the input. A dense layer will do the same operation, but acting only on a single grid cell of the input stack, acting at the same location upon all input feature maps within the stack at a time—using all values at that location (i.e., the 1 × 1 subsection) and transforming them according to assigned weights to pass forward a new set of channels to a single grid cell on the output stack. Each layer, either kernel or dense, may expand or contract the number of channels it passes forward. A kernel layer may also change the width and height dimensions of the subsection it passes forwards.We require an output that corresponds to a map of fire-occurrences; therefore our model needs to perform a series of transforms that preserve the width and height of the input, but reduce it to a single channel. The single channel in the output then represents the probability of each cell being classified as fire or not-fire (0–1).Our CNN model is comprised of 5 kernel layers (K1–K5 in Fig. 5), each acts on a 3 × 3 subsection and preserves width and height, passing forwards a transformed 3 × 3 section. Kernel K1 takes an input of 31 channels (predictor variables) but passes forward 128 channels to form the transformation T1 (Fig. 6). Kernels K2–K4 take inputs of 128 channels and pass forward 128 channels (T2–T4). Kernel K5 takes an input of 128 channels but passes forward 1 channel—the output. After each kernel applies its weights, there is an activation function applied before the values are passed on, which modify the answer to fit the necessary criteria to be a valid input to the next process. Kernels K1–K4 have a rectified linear (relu) activation function, which returns the input value if positive, and 0 if negative. Kernel K5 has a sigmoid activation function, that transforms the input values to between 0 and 1 such that negative values are transformed to 0.5.Fig. 6: Model structural diagram.Model structural diagram showing the input, 3 × 3 kernel layers (K1–K5), each transformation passed forwards (T1–T4) and the output, with all dimensions labelled.Full size imageModel training and validationWe used a stochastic gradient descent optimising function called Adam41 combined with a binary cross-entropy loss function to train the model against our fire-hotspot dataset iterated over 20 epochs. We split the data 70/30, using 70% as training data and 30% as validation data, recording accuracy, precision, and recall as the performance metrics, as well as the loss function itself.After model training, we applied the model to each yearly raster stack and compared the output against the fire-hotspot data for further model validation. Before validating the model outputs, we applied a simple 3 × 3 moving average window as a smoothing function to reduce the edge effects of tiling that are a by-product of having to split the study area into smaller tiles (32 × 32) for passing to the model. For this yearly validation, we again used the metrics accuracy, precision, and recall, such that:$${{{{{rm{Accuracy}}}}}}=100({{{{{rm{TP}}}}}}+{{{{{rm{TN}}}}}})/({{{{{rm{TP}}}}}}+{{{{{rm{TN}}}}}}+{{{{{rm{FP}}}}}}+{{{{{rm{FN}}}}}})$$$${{{{{rm{Precision}}}}}}=100({{{{{rm{TP}}}}}})/({{{{{rm{TP}}}}}}+{{{{{rm{FP}}}}}})$$$${{{{{rm{Recall}}}}}}=100({{{{{rm{TP}}}}}})/({{{{{rm{TP}}}}}}+{{{{{rm{FN}}}}}})$$where TP is true positive, TN is true negative, FP is false positive, and FN is false negative. These comparisons were made on a raster cell to raster cell basis after designating a 500 m buffer around each fire hotspot observation (MODIS and VIIRS data) and converting the buffers to a raster image of the same resolution and extent as the model prediction.ScenariosAfter validating the model performance, we built future scenarios to investigate the impact on fire occurrence of managing key anthropogenic features of the landscape: canals and land cover (Table 2).Table 2 Future scenario types and descriptions.Full size tableStudies have shown that unmanaged areas of heavily degraded or cleared swamp-forest are most susceptible to fires16,17,25,26,33,42. Therefore, we have built scenarios that investigate the possible impact of managing these areas by altering the model inputs to re-assign the land-cover designations ‘Swamp shrubland’ and ‘Scrubland’, as well as other land designation alterations. The first such restoration scenario investigates the impact of reforesting these areas by re-assigning the designations to ‘Swamp forest’. The second such scenario investigates the impact of converting these unmanaged areas to plantations by re-assigning the designations to ‘Plantation’. We also built two further land cover scenarios to investigate the impact of continued deforestation in the region by re-assigning the ‘Swamp forest’ designation to ‘Swamp shrubland’ and ‘Plantation’.We then built a scenario to investigate the impact of canal blocking on fire occurrence, modifying the proximity to canals model input by reducing the number of canals included in our proximity analysis to just two major canals, one that runs north-south, and one that runs west-east (Fig. 1). These canals could not practically be blocked due to their size and importance as navigation conduits.The final scenario simulates the combined impact of both re-foresting unmanaged degraded and cleared forest areas and the blocking of canals simultaneously.To evaluate the impact of each scenario on fire occurrences, we calculated the ratio of model predictions >0.5 probability (i.e., that a fire would occur in that raster cell) for each year for each scenario against the same year for the baseline scenario.Model use as a predictive toolTo evaluate the model’s potential to predict future fire distribution across the wider ex-Mega Rice Project area, we trained a second version of the model following the same methodology outlined above, but included only data from 2002 to 2018 in the training and test data passed to the model fitting algorithm. We then applied the model to the predictor variables corresponding to 2019 and compared model outputs to the observations of fire-occurrences by again looking at the metrics accuracy, precision, and recall. We also present a visual comparison of the outputs from the full model (2019 included in training data), the predictive model (2019 not included), and the observation data (MODIS and VIIRS hotspots). More

Taco-Taype, N. & Zúñiga-Dávila, D. Efecto de la inoculación de plantas de Tarwi con cepas de Bradyrhizobium spp. aisladas de un lupino silvestre, en condiciones de invernadero. Revista peruana de biología. 27, 35–42 (2022).Article 

Google Scholar 
Atchison, G. W. et al. Lost crops of the Incas: Origins of domestication of the Andean pulse crop tarwi Lupinus mutabilis. Am. J. Bot. 103, 1592–1606 (2016).CAS 
Article 

Google Scholar 
Peru Origins. Tarwi (Lupinus Mutabilis). https://peruorigins.com/tarwi/ (2022).Guilengue, N., Alves, S., Talhinhas, P. & Neves-Martins, J. Genetic and genomic diversity in a tarwi (Lupinus mutabilis Sweet) germplasm collection and adaptability to Mediterranean climate conditions. Agronomy 10, 21 (2020).Article 

Google Scholar 
Repo-Carrasco-Valencia, R., Basilio-Atencio, J., Luna-Mercado, G. I., Pilco-Quesada, S. & VidaurreRuiz, J. Andean ancient grains: Nutritional value and novel uses. Biol. Life Sci. Forum. https://doi.org/10.3390/blsf2021008015 (2022).Article 

Google Scholar 
Gulisano, A., Alves, S., Martins, J. N. & Trindade, L. M. Genetics and breeding of Lupinus mutabilis: An emerging protein crop. Front. Plant Sci. https://doi.org/10.3389/fpls.2019.01385 (2019).Article 

Google Scholar 
Chen, Y., She, Y., Zhang, R., Wang, J. & Zhang, X. Use of starch-based fat replacers in foods as a strategy to reduce dietary intake of fat and risk of metabolic diseases. Food Sci. Nutr. 8, 16–22 (2020).CAS 
Article 

Google Scholar 
Frick, K. M., Kamphuis, L. G., Siddique, K. H. M., Singh, K. B. & Foley, R. C. Quinolizidine alkaloid biosynthesis in lupins and prospects for grain quality improvement. Front. Plant Sci. https://doi.org/10.3389/fpls.2017.00087 (2017).Article 

Google Scholar 
Chirinos-Arias, M. C. Andean Lupin (Lupinus mutabilis Sweet) a plant with nutraceutical and medicinal potential. Revista Bio. Ciencias. 3, 163–172 (2015).
Google Scholar 
Wink, M. Chemical defense of leguminosae. Are quinolizidine alkaloids part of the antimicrobial defense system of lupins?. Zeitschrift für Naturforschung C. 39, 548–552 (1984).Article 

Google Scholar 
Hidalgo, M. et al. Evaluation of in vitro suceptibility to spartein in four strain of Mycobacterium tuberculosis. Rev. Peru Med Exp Salud Publica. 39, 77–82 (2022).Article 

Google Scholar 
Muñoz, E. B., Luna-Vital, D. A., Fornasini, M., Baldeón, M. E. & Gonzalez de Mejia, E. Gamma-conglutin peptides from Andean lupin legume (Lupinus mutabilis Sweet) enhanced glucose uptake and reduced gluconeogenesis in vitro. J. Funct. 45, 339–347 (2018).Article 

Google Scholar 
Bryant, L., Rangan, A. & Grafenauer, S. Lupins and health outcomes: A systematic literature review. Nutrients 14, 327 (2022).CAS 
Article 

Google Scholar 
Jacobsen, S. & Mujica, A. Geographical distribution of the Andean lupin (Lupinus mutabilis Sweet). Plant Genet. Resour. Newslett. 155, 1–8 (2008).
Google Scholar 
Antunez de Mayolo, S. Nutricion en el antiguo Peru. Banco Central de la Republica. Lima, Peru. 127 (1981).FAO. Perfiles nutricionales por paises: Peru. (ed. FAO) 36 p. (2000).UNICEF. Estado Mundial de la Infancia 2019 incluye a Perú entre las experiencias exitosas de lucha contra la desnutrición crónica infantile. https://www.unicef.org/peru/nota-de-prensa/estado-mundial-infancia-nutricion-alimentos-derechos-peru-experiencias-exitosas-desnutricion-cronica-infantil-reporte (2022).MINSA (Ministry of health – Peru). Situacion actual de la anemia. https://anemia.ins.gob.pe/situacion-actual-de-la-anemia-c1 (2022).WHO. Anemia. https://www.who.int/es/health-topics/anaemia#tab=tab_1 (2022).Galani, Y. J. H., Orfila, C. & Gong, Y. Y. A review of micronutrient deficiencies and analysis of maize contribution to nutrient requirements of women and children in Eastern and Southern Africa. Crit. Rev. Food Sci. Nutr. 62, 1568–1591 (2022).CAS 
Article 

Google Scholar 
White, P. J. & Martin, R. B. Biofortifying crops with essential mineral elements. Trends Plant Sci. 10, 586–593 (2005).Article 

Google Scholar 
White, P. J. & Martin, R. B. Biofortification of crops with seven mineral elements often lacking in human diets-iron, zinc, copper, calcium, magnesium, selenium and iodine. New Phytol. 182, 49–84 (2009).CAS 
Article 

Google Scholar 
Waters, B. M. & Sankaran, R. P. Moving micronutrients from the soil to the seeds: genes and physiological processes from a biofortification perspective. Plant Sciences. 180, 562–574 (2011).CAS 
Article 

Google Scholar 
Brooker, R. W. et al. Improving intercropping: A synthesis of research in agronomy, plant physiology and ecology. New Phytol. 206, 107–117 (2015).Article 

Google Scholar 
Ducsay, L. et al. Possibility of selenium biofortification of winter wheat grain. Plant Soil Environ. 62, 379–383 (2016).CAS 
Article 

Google Scholar 
Kumar, S. & Pandey, G. Biofortification of pulses and legumes to enhance nutrition. Heliyon. https://doi.org/10.1016/j.heliyon.2020.e03682 (2020).Article 

Google Scholar 
Diehn, T. A. et al. Boron demanding tissues of Brassica napus express specific sets of functional Nodulin26-like Intrinsic Proteins and BOR 1 transporters. Plant J. 100, 68–82 (2019).CAS 
Article 

Google Scholar 
Jayalakshmi, V. A., Reddy, T. & Nagamadhuri, K. V. Genetic diversity and variability for protein and micro nutrients in advance breeding lines and chickpea varieties grown in Andhra Pradesh.”. Legume Res. Int. J. 42, 768–772 (2019).
Google Scholar 
Bouis, H. & Saltzman, A. Improving nutrition through biofortification: A review of evidence from HarvestPlus, 2003 through 2016. Glob Food Sec. 12, 49–58 (2017).Article 

Google Scholar 
Sanca, D. Composición nutricional de diez genotipos de lupino (L. mutabilis y L. albus) desamargados por proceso acuoso. Thesis. Universidad Nacional Agraria La Molina. (2015).Rodríguez, A. Evaluación “in vitro” de la actividad antibacteriana de los alcaloides del agua de desamargado del chocho (Lupinus mutrabilis Sweet). Thesis. Escuela Superior Politécnica de Chimborazo, Ecuador (2009).Villacres, E. et al. Germination, an effective process to in-crease the nutritional value and reduce non-nutritive factors of lupine grain (Lupinus mutabilis Sweet). Int. J. Food Sci. Nutr. Eng. 5, 163–168 (2015).
Google Scholar 
Villacres, E., Rubio, A., Egas, L., Segovia, G. Usos alternativos del chocho: Chocho (Lupinus mutabilis Sweet) alimento andino redescubierto. IOP publishing: repositorio. https://repositorio.iniap.gob.ec/handle/41000/298 (2006).Ortega-David, E. A., Rodríguez, A. D. & Burbano, A. Z. Caracterización de semillas de lupino (Lupinus mutabilis) sembrado en los Andes de Colombia. Acta Agronómica. 59, 111–118 (2010).
Google Scholar 
White, P. J. & Broadley, M. R. Physiological limits to zinc biofortification of edible crops. Front Plant Sci. 80, 1–11 (2011).
Google Scholar 
Zhao, F., Su, Y. H., Dunham, S. J. & Rakszegiet, M. Variation in mineral micronutrient concentrations in grain of wheat lines of diverse origin. J. Cereal Sci. 49, 290–295 (2009).CAS 
Article 

Google Scholar 
Uauy, C., Distelfeld, A., Fahima, T., Blechl, A. & Dubcovsky, J. A NAC Gene regulating senescence improves grain protein, zinc, and iron content in wheat. Science 24, 1298–1301 (2006).ADS 
Article 

Google Scholar 
Shorrocks, V. M. The occurrence and correction of boron deficiency. Plant Soil 193, 121–148 (1997).CAS 
Article 

Google Scholar 
D’Imperio, M. et al. Boron biofortification of Portulaca oleracea L. through soilless cultivation for a new tailored crop. Agronomy. 10, 999–1013 (2020).Article 

Google Scholar 
Boyacioglu, O., Orenay-Boyacioglu, S., Yildirim, H. & Korkmaz, M. Boron intake, osteocalcin polymorphism and serum level in postmenopausal osteoporosis. J. Trace Elem. Med. Biol. 48, 52–56 (2018).CAS 
Article 

Google Scholar 
Oliveira Araújo, E., Ferreira Dos Santos, E. & Camacho Oliveira, M. A. Boron-zinc interaction in the absorption of micronutrients by cotton. Agronomía Colombiana. 36, 51–57 (2018).Article 

Google Scholar 
Squitti, R., Siotto, M. & Polimanti, R. Low-copper diet as a preventive strategy for Alzheimer’s disease. Neurobiol. Aging 2, 40–50 (2014).Article 

Google Scholar 
Schilsky, M.L. Management of Wilson Disease (A Pocket Guide), 1st ed.; Publisher: Humana Press, Farmington, CT, USA. 154–196 (2018).Martins, A. C. et al. Manganese in the diet: Bioaccessibility, adequate intake, and neurotoxicological effects. J. Agric. Food Chem. 46, 12893–12903 (2020).Article 

Google Scholar 
Falah, S. A. & Saja, N. M. Essential trace elements and their vital roles in human body. Indian J. Adv. Chem. Sci. 3, 127–136 (2017).
Google Scholar 
National institutes of health. Manganese. Fact Sheet for Health Professionals. IOP Publishing ods.od.nih.gov. https://ods.od.nih.gov/factsheets/Manganese-HealthProfessional/. (2021).Savadi, S. Molecular regulation of seed development and strategies for engineering seed size in crop plants. Plant Growth Regul. 84, 401–422 (2018).CAS 
Article 

Google Scholar 
Ge, L. et al. (2016) Increasing seed size and quality by manipulating BIG SEEDS1 in legume species. Proc Natl Acad Sci. 113, 12414–12419 (2016).CAS 
Article 

Google Scholar 
Zou, L. Effects of gradual and sudden heat stress on seed quality of Andean lupin, Lupinus mutabilis. Thesis. University of Helsinki. https://helda.helsinki.fi/handle/10138/16501 (2009).Buircell, B.J., Cowling, A.W. Genetic Resources in Lupins (eds. Gladstones, J.S., Atkins, C.A., Hamblin, J.) (United Kingdom: CAB International, 1998).Aguilar-Angulo, L. A. Evaluación del rendimiento de grano y capacidad simbiótica de once accesiones de tarwi (Lupinus mutabilis Sweet), bajo condiciones de Otuzco-La Libertad (Universidad Nacional Agraria La Molina, 2015).
Google Scholar 
De La Cruz, N. Caracterización fenotípica y de rendimiento preliminar de ecotipos de tarwi (Lupinus mutabilis sweet), bajo condiciones del Callejón de Huaylas – Ancash (Universidad Nacional Agraria la Molina, 2018).
Google Scholar 
Huisa, J. Evaluación del comportamiento agronómico de catorce accesiones del ensayo nacional de tarwi (Lupinus mutabilis sweet.) en el CIP Camacani Puno – Perú”. Thesis. Universidad Nacional Agraria la Moina (2018).Cayo, B. Evaluación del comportamiento agronómico de ocho genotipos selectos de tarwi (Lupinus mutabilis sweet) bajo condiciones del CIP. CAMACANI – UNA – PUNO. Thesis. Universidad Nacional del Altiplano (2020).Buircell, B.J., Cowling, A.W. Lupin. Lupinus spp. Promoting the conservation and use of underutilized and ne-glected crops (eds. Gladstones, J.S., Atkins, C.A., Hamblin, J.) (United Kingdom: CAB International, 1998).Plata, J. Comportamiento Agronómico de dos Variedades de tarwi (Lupinus mutabilis Sweet), bajo tres densidades de siembra en la comunidad Marka Hilata Carabuco (Universidad San Andres, 2016).
Google Scholar 
Mendoza, C. Rendimiento de ecotipos regionales y variedades de tarwi (Lupínus mutabilis Sweet.) en el valle del Mantaro, Jauja, Junín. Thesis. Universidad Nacional Agraria la Moina (2020).Aguilar, S. Sistemas de producción de Lupinus mutabilis Sweet ‘chocho’ en terrazas y laderas con fertilización fosfatada en Cajamarca. Dissertation. La Molina National Agrarian University (2011).Aquino, S. Sustentabilidad del cultivo de tarwi (Lupinus mutabilis sweet) en la zona altoandina del Valle del Mantaro (Universidad Nacional Agraria la Molina, 2018).
Google Scholar 
Barda, M. S., Chatzigeorgiou, T., Papadopoulos, G. K. & Bebeli, P. J. Agro-morphological evaluation of Lupinus mutabilis in two locations in greece and association with insect pollinators. Agriculture https://doi.org/10.3390/agriculture11030236 (2021).Article 

Google Scholar 
Herniter, I. A., Jia, Z. & Kusi, F. Market preferences for cowpea (Vigna unguiculata [L.] Walp) dry grain in Ghana. African J Ag Res. 14, 928–934 (2019).Article 

Google Scholar 
Dordas, C. Foliar boron application affects lint and seed yield and improves seed quality of cotton grown on calcareous soils. Nutr. Cycl. Agroecosyst. 76, 19–28 (2006).CAS 
Article 

Google Scholar 
Kristek, S. et al. Effect of various rates of boron on yield and quality of high-grade sugar beet varieties. Listy Cukrovarnické a Řepařské. 4, 146–150 (2018).
Google Scholar 
Thomas, C. L. et al. Root morphology and seed and leaf ionomic traits in a Brassica napus L. diversity panel show wide phenotypic variation and are characteristic of crop habit. BMC Plant Biol. 16, 214–232 (2016).CAS 
Article 

Google Scholar 
Dursun, A. et al. Effects of boron fertilizer on tomato, pepper and cucumber yields and chemical composition. Commun Soil Sci Plant Anal. 1, 1576–1593 (2010).Article 

Google Scholar 
Sotiropoulos, T. E., Therios, T. N., Dimassi, K. N., Bosabalidis, A. & Kofidis, G. Nutritional status, growth, CO2 assimilation, and leaf anatomical responses in two kiwifruit species under boron toxicity. J Plant Nutr. 25, 1249–1261 (2002).CAS 
Article 

Google Scholar 
Muccifora, S. & Bellani, L. Effects of copper on germination and reserve mobilization in Vicia sativa L. seeds. Environ. Pollut. 179, 68–74 (2013).CAS 
Article 

Google Scholar 
Kobraee, S. Effect of foliar fertilization with zinc and manganese sulfate on yield, dry matter accumulation, and zinc and manganese contents in leaf and seed of chickpea (Cicer arietinum). J. Appl. Biol. Biotechnol. 7, 20–28 (2019).CAS 

Google Scholar 
IBPGR (1981) Lupin descriptors. https://www.bioversityinternational.org/fileadmin/bioversity/publications/Web_version/103/ (1981).Zasoski, R. J. & Burau, R. G. A rapid nitric-perchloric acid digestion method for multi-element tissue analysis. Commun. Soil Sci. Plant Anal. 8, 425–436 (1997).Article 

Google Scholar 
Pereira, T., Coelho, C. M. M., Bogo, A., Guidolin, A. F. & Miquelluti, D. J. Diversity in common bean landraces from south Brazil. Acta Bot. Croat. 1, 79–92 (2009).
Google Scholar 
Pujar, M., Govindaraj, M., Gangaprasad, S., Kanatti, A. & Shivade, H. Genetic variation and diversity for grain iron, zinc, protein and agronomic traits in advanced breeding lines of pearl millet [Pennisetum glaucum (L.) R Br] for biofortification breeding. Genet. Resour. Crop Evol. 67, 2009–2022 (2020).CAS 
Article 

Google Scholar 
Lira, J. P. E. et al. Safflower genetic diversity based on agronomic characteristics in Mato Grosso state, Brazil, for a crop improvement program. Genet. Mol. Res. 1, 1–12 (2021).
Google Scholar 
de Sá, S. F. et al. Genetic diversity via REML-BLUP of ex situ conserved macauba [Acrocomia aculeata (Jacq.) Lodd. ex Mart.] ecotypes. Genet. Resour. Crop Evol. 68, 3193–3204 (2021).Article 

Google Scholar 
Kuru, R., Yilmaz, S., Tasli, P. N., Yarat, A. & Sahin, F. Boron content of some foods consumed in Istanbul, Turkey. Biol. Trace Elem. Res. 187, 1–8 (2019).CAS 
Article 

Google Scholar 
Shokunbi, O., Adepoju, O., Mojapelo, P., Ramaite, I. & Akinyele, I. Copper, manganese, iron and zinc contents of Nigerian foods and estimates of adult dietary intakes. J. Food Compos. Anal. 82, 103–245 (2019).Article 

Google Scholar 
Norwegian scientific committee for food and environment. Assessment of dietary intake of manganese in rela-tion to tolerable upper intake. IOP Publishing wkm. www.vkm.no. (2018).Gil, V., Guzmán, L. & Quintero, E. Caracterización de la variabilidad morfológica de un “genotipo local” de maíz y dos de sus selecciones. Centro Agrícola. 4, 79–83 (2004).
Google Scholar  More

In early 2020, my colleagues and I realized that animal-tracking data collected before, during and after the pandemic lockdowns could provide invaluable insights into human–wildlife interactions and conservation benefits on a global scale. We launched a research consortium — the COVID-19 Bio-Logging Initiative — to investigate how animals behaved while much of the world’s human population sheltered at home.But we had no way to establish how many, and which, animals were wearing tags. Miniature tracking devices are routinely attached to a vast range of species — from songbirds to whales — to collect detailed data on their movements, behaviour and physiology. Yet, of the thousands of ‘bio-loggers’ deployed every year, many generate data sets that remain effectively undiscoverable — they are saved on personal hard drives or institutional servers, inaccessible to the wider community. This problem can be solved by setting up a global registry for all tags on wild animals.Although individual tracking studies make important contributions to our understanding of the ecological needs of animal species, pooling data (across taxa, longer time periods or multiple locations) can reveal general patterns, aiding the design of particularly effective conservation strategies. For example, integrating the tracks of 4,060 animals across 17 marine species (including albatrosses, penguins, seals and whales) has helped to identify conservation priority areas in the Southern Ocean (M. A. Hindell et al. Nature 580, 87–92; 2020).In an ideal world, all animal-tracking data would be archived — with either open or restricted access — in public repositories, such as Movebank. Excellent progress has been made towards this goal, but universal uptake is hindered by time constraints, governmental or institutional restrictions and concerns over inappropriate data use.To encourage as many data owners as possible to join the COVID-19 Bio-Logging Initiative, we launched a recruitment campaign through Movebank, social media, mailing lists, newsletters, personal contacts and a published call to action (C. Rutz et al. Nature Ecol. Evol. 4, 1156–1159; 2020). Our consortium has grown to more than 600 international collaborators, accumulating a staggering one billion location records for some 200 animal species. Despite this impressive community response, we know that this is only the tip of the iceberg.The global tag registry that I suggest would contain metadata for tags (including tag type and settings, information on the animal, and date and location of deployment), as well as researchers’ contact details — but not the actual tracking data. This decoupling of information would unlock the field’s full conservation potential in the short term and would build the trust required to allow raw data to be archived routinely in public repositories in the longer term. Over time, the tag registry is likely to evolve naturally into a ‘meta-repository’, linking to raw data sets hosted across a multitude of repositories.The registry would enable researchers to check data availability at the push of a button — for example, for a particular taxonomic group, such as terrestrial carnivores, or a specific region, such as the Pacific Ocean — and to get in touch with the relevant data owners. Registry management must comply with international best practices, so robust processes would need to be set up to vet queries, pass on collaboration proposals to data owners and minimize overlap between studies.For the registry to fulfil its intended purpose, it must be used by the entire animal-tracking community. How can this be achieved? I see an opportunity to integrate tag registration into existing ethical-review processes. Governmental authorities, research institutions, funders, publishers and fieldworkers agree that permits must be in place before animals can be tagged. Building on this international consensus, ethical review boards could make tag registration a condition of study approval.To complement this bottom-up approach, well established initiatives — such as those associated with the United Nations Environment Programme or the International Union for Conservation of Nature — could help to build an international policy mandate and provide independent oversight. The International Bio-Logging Society, which has been working to unite animal-tracking efforts on land and at sea, could provide crucial support.This vision is no doubt ambitious, but it is achievable. Every civil aircraft on the planet must be registered — so I am convinced that, with effective coordination, we can accomplish the same for tagged animals. Furthermore, the basic principle of hosting metadata, but not raw data, is being used productively by other databases, such as AviSample — a registry for biological samples collected from wild birds.Many researchers, myself included, feel a moral obligation to the animals carrying our tags. A global tag registry would help to realize the full conservation potential of all tracking data, minimize duplication of tagging efforts and facilitate sharing of welfare-related expertise. The conservation cost of missing data in large-scale collaborative projects cannot be easily measured, but is probably substantial. We simply cannot afford this, and must ensure that all animal-tracking data are immediately discoverable.

Competing Interests
This article is a contribution of the COVID-19 Bio-Logging Initiative, which is funded in part by the Gordon and Betty Moore Foundation (GBMF9881) and the National Geographic Society (NGS-82515R-20) (both grants to C.R.), and endorsed by the United Nations Decade of Ocean Science for Sustainable Development. More

Chaigneau, A., Gizolme, A. & Grados, C. Mesoscale eddies off Peru in altimeter records: identification algorithms and eddy spatio-temporal patterns. Prog. Oceanogr. 79, 106–119 (2008).ADS 
Article 

Google Scholar 
McGillicuddy, D. J. Jr et al. Influence of mesoscale eddies on new production in the Sargasso Sea. Nature 394, 263–266 (1998).ADS 
CAS 
Article 

Google Scholar 
Dufois, F. et al. Anticyclonic eddies are more productive than cyclonic eddies in subtropical gyres because of winter mixing. Sci. Adv. 2, 1–7 (2016).Article 

Google Scholar 
Godø, O. R. et al. Mesoscale eddies are oases for higher trophic marine life. PLoS ONE 7, e30161 (2012).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Chelton, D. B., Gaube, P., Schlax, M. G., Early, J. J. & Samelson, R. M. The influence of nonlinear mesoscale eddies on near-surface oceanic chlorophyll. Science 334, 328–333 (2011).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Sarmiento, J. L. et al. Response of ocean ecosystems to climate warming. Global Biogeochem. Cycles 18, GB3003 (2004).ADS 
Article 
CAS 

Google Scholar 
Bell, J. D. et al. Diversifying the use of tuna to improve food security and public health in Pacific Island countries and territories. Mar. Policy 51, 584–591 (2015).Article 

Google Scholar 
Della Penna, A. & Gaube, P. Mesoscale eddies structure mesopelagic communities. Front. Mar. Sci. 7, 454 (2020).ADS 
Article 

Google Scholar 
Braun, C. D. et al. The functional and ecological significance of deep diving by large marine predators. Ann. Rev. Mar. Sci. 14, 129–159 (2022).PubMed 
Article 

Google Scholar 
McGillicuddy, D. J. Jr Mechanisms of physical-biological-biogeochemical interaction at the oceanic mesoscale. Ann. Rev. Mar. Sci. 8, 125–159 (2016).PubMed 
Article 

Google Scholar 
Fennell, S. & Rose, G. Oceanographic influences on deep scattering layers across the North Atlantic. Deep-Sea Res. Part I Oceanogr. Res. Pap. 105, 132–141 (2015).ADS 
Article 

Google Scholar 
Duffy, L. M. et al. Global trophic ecology of yellowfin, bigeye, and albacore tunas: understanding predation on micronekton communities at ocean-basin scales. Deep-Sea Res. Part II Topical Stud. Oceanogr. 140, 55–73 (2017).ADS 
Article 

Google Scholar 
Gaube, P. et al. Mesoscale eddies influence the movements of mature female white sharks in the Gulf Stream and Sargasso Sea. Sci. Rep. 8, 7363 (2018).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Braun, C. D., Gaube, P., Sinclair-Taylor, T. H., Skomal, G. B. & Thorrold, S. R. Mesoscale eddies release pelagic sharks from thermal constraints to foraging in the ocean twilight zone. Proc. Natl Acad. Sci. USA 116, 17187–17192 (2019).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Doyle, T. K. et al. Leatherback turtles satellite-tagged in European waters. Endanger. Species Res. 4, 23–31 (2008).Article 

Google Scholar 
Pauly, D. & Christensen, V. Primary production required to sustain global fisheries. Nature 374, 255–257 (1995).ADS 
CAS 
Article 

Google Scholar 
Lynham, J., Nikolaev, A., Raynor, J., Vilela, T. & Villaseñor-Derbez, J. C. Impact of two of the world’s largest protected areas on longline fishery catch rates. Nat. Commun. 11, 979 (2020).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Polovina, J. J., Abecassis, M., Howell, E. A. & Woodworth, P. Increases in the relative abundance of mid-trophic level fishes concurrent with declines in apex predators in the subtropical North Pacific, 1996-2006. Fish. Bull. 107, 523–531 (2009).
Google Scholar 
Royer, T. C. Ocean eddies generated by seamounts in the North Pacific. Science 199, 1063–1064 (1978).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Liu, Y. et al. Eddy analysis in the subtropical zonal band of the North Pacific Ocean. Deep-Sea Res. Part I Oceanogr. Res. Pap. 68, 54–67 (2012).ADS 
Article 

Google Scholar 
Bernstein, R. L. & White, W. B. Time and length scales of baroclinic eddies in the central North Pacific Ocean. J. Phys. Oceanogr. 4, 613–624 (1974).ADS 
Article 

Google Scholar 
Maunder, M. N. & Punt, A. E. Standardizing catch and effort data: a review of recent approaches. Fish. Res. 70, 141–159 (2004).Article 

Google Scholar 
Woodworth, P. A. et al. Eddies as offshore foraging grounds for melon-headed whales (Peponocephala electra). Mar. Mammal Sci. 28, 638–647 (2012).Article 

Google Scholar 
Gaube, P. et al. The use of mesoscale eddies by juvenile loggerhead sea turtles (Caretta caretta) in the southwestern Atlantic. PLoS ONE 12, e0172839 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Chambault, P. et al. Swirling in the ocean: immature loggerhead turtles seasonally target old anticyclonic eddies at the fringe of the North Atlantic Gyre. Prog. Oceanogr. 175, 345–358 (2019).ADS 
Article 

Google Scholar 
Gaube, P., McGillicuddy Jr, D., Chelton, D., Behrenfeld, M. & Strutton, P. Regional variations in the influence of mesoscale eddies on near-surface chlorophyll. J. Geophys. Res. Oceans 119, 8195–8220 (2014).Waga, H., Kirawake, T. & Ueno, H. Impacts of mesoscale eddies on phytoplankton size structure. Geophys. Res. Lett. 46, 13191–13198 (2019).ADS 
Article 

Google Scholar 
Irigoien, X. et al. Large mesopelagic fishes biomass and trophic efficiency in the open ocean. Nat. Commun. 5, 3271 (2014).ADS 
PubMed 
Article 
CAS 

Google Scholar 
Chen, Y.-lL. et al. Biologically active warm-core anticyclonic eddies in the marginal seas of the western Pacific Ocean. Deep Sea Res. Part I 106, 68–84 (2015).CAS 
Article 

Google Scholar 
Harke, M. J. et al. Microbial community transcriptional patterns vary in response to mesoscale forcing in the North Pacific Subtropical Gyre. Environ. Microbiol. 23, 4807–4822 (2021).CAS 
PubMed 
Article 

Google Scholar 
Hawco, N. J. et al. Iron depletion in the deep chlorophyll maximum: mesoscale eddies as natural iron fertilization experiments. Global Biogeochem. Cycles 35, e2021GB007112 (2021).ADS 
CAS 
Article 

Google Scholar 
Klevjer, T. A. et al. Large scale patterns in vertical distribution and behaviour of mesopelagic scattering layers. Sci. Rep. 6, 19873 (2016).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Behrenfeld, M. J. et al. Global satellite-observed daily vertical migrations of ocean animals. Nature 576, 257–261 (2019).CAS 
PubMed 
Article 

Google Scholar 
Madigan, D. J. et al. Water column structure defines vertical habitat of twelve pelagic predators in the South Atlantic. ICES J. Mar. Sci. 78, 867–883 (2021).Article 

Google Scholar 
Arostegui, M., Gaube, P. & Braun, C. Movement ecology and stenothermy of satellite-tagged shortbill spearfish (Tetrapturus angustirostris). Fish. Res. 215, 21–26 (2019).Article 

Google Scholar 
Lehodey, P., Senina, I. & Murtugudde, R. A spatial ecosystem and populations dynamics model (SEAPODYM)—modeling of tuna and tuna-like populations. Prog. Oceanogr. 78, 304–318 (2008).ADS 
Article 

Google Scholar 
Varghese, S. P., Somvanshi, V. S. & Dalvi, R. S. Diet composition, feeding niche partitioning and trophic organisation of large pelagic predatory fishes in the eastern Arabian Sea. Hydrobiologia 736, 99–114 (2014).CAS 
Article 

Google Scholar 
Ward, P. & Myers, R. A. Inferring the depth distribution of catchability for pelagic fishes and correcting for variations in the depth of longline fishing gear. Can. J. Fish. Aquat.Sci. 62, 1130–1142 (2005).Article 

Google Scholar 
Kai, E. T. et al. Top marine predators track Lagrangian coherent structures. Proc. Natl Acad. Sci. USA 106, 8245–8250 (2009).ADS 
CAS 
Article 

Google Scholar 
Lima, I. D., Olson, D. B. & Doney, S. C. Biological response to frontal dynamics and mesoscale variability in oligotrophic environments: biological production and community structure. J. Geophys. Res. Oceans 107, 25-1–25-21 (2002).Article 

Google Scholar 
Spall, S. A. & Richards, K. J. A numerical model of mesoscale frontal instabilities and plankton dynamics—I. model formulation and initial experiments. Deep-Sea Res. Part I Oceanogr. Res. Pap. 47, 1261–1301 (2000).ADS 
Article 

Google Scholar 
Siegelman, L., O’Toole, M., Flexas, M., Rivière, P. & Klein, P. Submesoscale ocean fronts act as biological hotspot for southern elephant seal. Sci. Rep. 9, 5588 (2019).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Lévy, M., Ferrari, R., Franks, P. J., Martin, A. P. & Rivière, P. Bringing physics to life at the submesoscale. Geophys. Res. Lett. https://doi.org/10.1029/2012GL052756 (2012).Article 

Google Scholar 
Guidi, L. et al. Does eddy-eddy interaction control surface phytoplankton distribution and carbon export in the North Pacific Subtropical Gyre? J. Geophys. Res. Biogeosciences https://doi.org/10.1029/2012JG001984 (2012).Article 

Google Scholar 
Chow, C. H., Cheah, W., Tai, J. H. & Liu, S. F. Anomalous wind triggered the largest phytoplankton bloom in the oligotrophic North Pacific Subtropical Gyre. Sci. Rep. 9, 15550 (2019).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Guo, M., Xiu, P., Chai, F. & Xue, H. Mesoscale and submesoscale contributions to high sea surface chlorophyll in subtropical gyres. Geophys. Res. Lett. 46, 13217–13226 (2019).ADS 
Article 

Google Scholar 
Klein, P. et al. Ocean-scale interactions from space. Earth Space Sci. 6, 795–817 (2019).ADS 
Article 

Google Scholar 
Martin, A. et al. The oceans’ twilight zone must be studied now, before it is too late. Nature 580, 26–28 (2020).ADS 
CAS 
PubMed 
Article 

Google Scholar 
St. John, M. A. et al. A dark hole in our understanding of marine ecosystems and their services: perspectives from the mesopelagic community. Front. Marine Sci. 3, 31 (2016).
Google Scholar 
Bigelow, K., Musyl, M. K., Poisson, F. & Kleiber, P. Pelagic longline gear depth and shoaling. Fish. Res. 77, 173–183 (2006).Article 

Google Scholar 
Brodziak, J. & Walsh, W. A. Model selection and multimodel inference for standardizing catch rates of bycatch species: a case study of oceanic whitetip shark in the Hawaii-based longline fishery. Can. J. Fish. Aquat.Sci. 70, 1723–1740 (2013).Article 

Google Scholar 
Woodworth-Jefcoats, P. A., Polovina, J. & Drazen, J. Synergy among oceanographic variability, fishery expansion, and longline catch composition in the central North Pacific Ocean. Fish. Bull. 116, 228–239 (2018).Article 

Google Scholar 
Boggs, C. H. Depth, capture time, and hooked longevity of longline-caught pelagic fish: timing bites of fish with chips. Fish. Bull. 90, 642–658 (1992).
Google Scholar 
Walsh, W. A. & Brodziak, J. Applications of Hawaii longline fishery observer and logbook data for stock assessment and fishery research. NOAA Tech. Memo. 57, 62 (2016).
Google Scholar 
Walsh, W. A. & Brodziak, J. Billfish CPUE standardization in the Hawaii longline fishery: model selection and multimodel inference. Fish. Res. 166, 151–162 (2015).Article 

Google Scholar 
Gilman, E., Chaloupka, M., Fitchett, M., Cantrell, D. L. & Merrifield, M. Ecological responses to blue water MPAs. PLoS ONE 15, e0235129 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Portner, E. J., Polovina, J. J. & Choy, C. A. Patterns in micronekton diversity across the North Pacific Subtropical Gyre observed from the diet of longnose lancetfish (Alepisaurus ferox). Deep-Sea Research Part I 125, 40–51 (2017).ADS 
Article 

Google Scholar 
Brooks, M. E. et al. glmmTMB balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling. R J. 9, 378–400 (2017).Article 

Google Scholar 
Hartig, F. DHARMa: Residual diagnostics for hierarchical (multi-level/mixed) regression models. R package version 0.3.3.0 http://florianhartig.github.io/DHARMa/ (2020).Jackson, C. H. Multi-state models for panel data: the msm package for R. J. Stat. Softw. https://doi.org/10.18637/jss.v038.i08 (2011).Article 

Google Scholar 
Bates, D. et al. lme4: Linear mixed-effects models using ’Eigen’ and S4. R package version 1.1-25 https://github.com/lme4/lme4/ (2020).Lenth, R. et al. emmeans: Estimated marginal means, aka least-squares mean. R package version 1.7.2 https://github.com/rvlenth/emmeans (2022).R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2020); http://www.r-project.org/ More

We used two agent-based simulation models to investigate the concepts of “cooperate for the spread” and “pay for the escape,” both were net logo models created by Dr. Susan Hanisch.Afterward, we modified the first model to represent the concept of sharing the dispersal costs. We used the second model without modifications. Instead, we assigned definite values of some parameters that highlight the pay for the escape strategy.First modelThe original model was entitled “Evolution and patchy resource”18. She first developed it for educational purposes. It illustrates the concepts of cooperator-cheater competition, natural selection, spatial structure mechanisms, multilevel selection, and founder effects.Changeable variables

Distance-resource-areas: the distance between the centers of the resource areas.

Size-resource areas: the size of resource areas as a radius in the number of patches.

Living costs: the costs that each agent has to deduct from energy per iteration for basic survival.

Mutation rate: The probability that offspring agents have different traits than their parents.

Evolution: the ability of agents to produce offspring.

Constant variables

The number of patches is 112 × 112 patches.

Carrying capacity per patch: Resource = 10, Agents = 1

The growth rate of the resource = 0.2

The resources on a patch regrow by a logistic growth function up to the carrying capacity: New resource level = current resource level + (Growth-Rate × current resource level) × (1 – (Current resource level/carrying capacity)).

The cost for producing offspring is ten subtracted units of energy.

The initial level of energy of agents is set at living costs.

Role of randomness

Agents are distributed randomly in resource areas at the beginning of a simulation.

Sustainable behavior is distributed randomly with a probability of percent sustainables among the initial agent population.

The order in which agents move and harvest within one iteration is random.

Agents move to a randomly selected patch if several patches fulfill the objectives.

The order in which agents produce offspring within one iteration is random.

Agents reproduce offspring with a probability of (0.0005 × Energy).

Agents place offspring on a randomly selected unoccupied neighboring patch.

Offspring mutate with a potential mutation rate.

Model processesIn each iteration, each agent moves around in random order. There are three likelihoods:

If there are no unoccupied patches in a two-patch radius, they stay on the current patch.

If there are unoccupied patches with resources amounting to more than living costs, the agents move to them.

If the resource amount is less than the living costs, the agents move randomly to other unoccupied patches.

The agents harvest the resources from separated patches to gain energy for metabolism and proliferation. If the energy level of any agent falls to zero, it dies. The cooperator type harvests half of the resource, while the greedy type consumes 99%.The living costs are deducted from the energy amount of the agent constantly everywhere all the time. This process occurs whether an agent moves within the patch, between the patches, or even not. Therefore, the model does not consider dispersal cost explicitly.If there is an unoccupied neighbor patch, the agent can reproduce with a probability of 0.0005 of his energy, place the offspring on the unoccupied neighbor patch, and then transfer ten units of the energy to his offspring.Resources regrow only on resource patches. When the resource amount is more than or equal to 0.1, then it regrows. When the resource is less than 0.1, its value is set to 0.1.Output diagrams and monitors

The average energy of agents: average energy levels of sustainable and greedy agents, resulting from resource harvest minus living costs and reproduction.

Trait frequencies: the relative frequencies of sustainable and greedy agents in the total population, resulting from mutations, different reproduction rates, and death.

Agent population: the absolute number of the total population size resulting from reproduction and death.

ModificationsIn the first modification, we added a different type of cost that agents only incur when they disperse from one patch to another (in-between the patches). It is the slider entitled “dispersal costs”.In the second modification, we added another sharing dispersal costs tool to reduce them by dividing their value by the number of included agents (flock-mates) in the identified range from the same type. It is the slider entitled “group-dispersal-range.” which is the flock mate’s areas as a radius in the number of patches. Therefore, changing the value of the group dispersal range will change the area around every agent. Accordingly, the number of its flock mates who share the dispersal costs also adjusts.The group dispersal range is not confined to greedy agents but applies to all agents. Therefore, it represents the case of the wild-type cooperators who can also cooperate for the spread. The group dispersal range also does not only target the agents in between patches. However, it counts the agents inside and outside the patches. For example, once an agent starts its dispersion with a determined range containing ten agents, four from another type, three non-dispersal agents from the same type that existed inside a patch, and three dispersal agents from the same type outside the patches. The dispersal costs for this agent will be divided by 6.Our assumption that non-dispersal agents at the pre-departure stage share dispersion costs with dispersal agents; seems justified because they reap mutual benefits by reducing kin competition inside patches if they promote the migrators. However, can agents remotely pay the dispersion costs? Yes. For instance, some bacterial species can trigger the migration of other species if located in their vicinity, even if the two bacterial colonies are separated by a barrier19,20 or if they are non-motile21. On the other hand, dispersion is an extended process with many factors, including escape from predators, suppression of host defense mechanisms, and production of biosurfactants to reduce surface tension to facilitate motility. Therefore, the agent’s contribution (inside/outside the patches) to support such factors is considered a shared dispersal cost.Finally, cheaters can arise within cooperator patches by mutation or immigration. Therefore, to investigate the efficacy of migration, the mutation rate value should be 0 to cancel its effect in the meta-population dynamics.Second modelThe model is entitled “Evolution, resources, monitoring, and punishment.”22 is a simulation of a population with four types of agents competing for the same resource. It demonstrates many concepts, such as kin selection, cooperation, selfishness, public good, monitoring, punishment, sharing the costs, positive/negative frequency-dependent selection, and multilevel selection. The four agent colors and types: (1) Red: greedy, non-punishing. (2) Orange: greedy, punishing. (3) Turquoise: sustainable, non-punishing. (4) Green: sustainable, punishing.Punishing agents can perceive other agents in their environment to some degree (perception accuracy) and react to their behavior. There are three kinds of punishment: Punishers can kill agents with greedy harvesting behavior, stop them from harvesting in the next iteration, or have them pay a penalty fee to their neighbors.Agents have a cost (energy) to pay for, both detection and punishment, so this behavior is altruistic. Punisher agents of one type share punishment costs equally.Changeable variables

Death rate: the probability that agents die independent of their energy level.

Carrying capacity: the maximum amount of resource units on a patch from 1 to 100.

Growth rate: the rate at which resources on patches regrow. The maximum sustainable yield is calculated based on the carrying capacity and growth rate.

Harvest-sustainable: the number of resource units harvested by sustainable agents.

Harvest-greedy: the number of resource units harvested by sustainable agents.

Perception accuracy: the probability with which punishing agents notice greedy agents.

Costs-perception: the costs in units of energy, punishing agents have to pay for perceiving other agents.

Costs-punishment: the costs as units of energy that punishing agents have to pay in each iteration to punish other agents. All punishing agents of an agent divide the costs of punishment.

Punishment: the kinds of punishing behavior that punishing agents perform.

Fine: if the kind of punishment is “pay fine”, the fine in energy units that punished agents have to pay (shared between all their neighbors).

Living costs and mutation rate: see the first model.

Constant variables

The number of patches: There are 60 × 60 patches in the world.

The initial energy level of agents is set at living costs + 1.

The initial number of resource units on a patch is set to the carrying capacity.

The resources on a patch regrow: see the first model.

Role of randomness* In addition to items in the first model.

Agents take on their traits (harvest preference and ability to notice and punish) randomly based on the probability of percent-sustainable and percent-punishers.

The order in which punishing agents notice greedy agents within one iteration is random.

Greedy agents are noticed by punishing agents with a probability of perception accuracy.

The order in which detected greedy agents are punished within one iteration is random.

Agents produce offspring with a probability of (0.001 × Energy).

Agents die with a probability of (death-rate).

Model processesIn each iteration, each agent attempts to harvest resources from the patches it is on and the eight neighboring patches until the harvest preference level is reached, except for the punished agent with the sanction (suspend harvest once), its harvest amount = 0 in the current iteration. If the amount of resources available is lower than the amount that the unpunished agent attempts to harvest. Then, the agent moves to a neighboring unoccupied patch with the most resources after losing one energy unit as a move cost.Punishers pay the costs of perceiving the greedy agents. The greedy neighbors have been noticed with the probability of perception accuracy. The agent lost an amount of energy as living costs. The agent dies with the likelihood of death rate or if the energy level falls to zero.If there is an unoccupied neighbor patch, the agent can reproduce with a probability of 0.001 of its energy, place the offspring on the unoccupied neighbor patch, and then transfer half of its energy to its offspring that mutate according to the probability of the mutation rate.Resources regrow on all patches. When the resource amount is more than or equal to 0.1, then it regrows. When the resource is less than 0.1, its value is set to 0.1.Output diagrams and monitors

Populations (% of carrying capacity): the state of the resource and the agent population in the world as a percentage of total carrying capacity resulting from resource harvesting behavior and resource regrowth, agent reproduction, and death.

Average harvest per iteration: the average harvested amounts of agents per iteration by trait, resulting from harvested resource units, minus costs for monitoring and punishing (for punishing agents), minus fines (for punished agents in case of punishment “Pay fine”)

The average energy of agents and trait frequencies: see the first model.

How does the model represent a conditional defector strategy?The model aims to highlight the role of kin selection and punishment mechanisms in supporting cooperation evolution against cheats. We did not need to modify the model but just thought about what the conditional defector should do to upside down the game. The answer was to pay for the escape.For instance, if the standard Harvest-greedy of a cheater (greedy, non-punishing) was 13 and the Perception-accuracy of its actual punishers was 75%. Now suppose this cheater faces troubles, and it cannot dominate. However, if it gives up some of its profit to become 12, to escape punishment, and to reduce the perception accuracy to 60%, it could dominate and take over the population.The conditional cheater can pay something and reduce its profit to escape punishment by reducing perception accuracy if there is a positive correlation between these two variables. Therefore, this model is appropriate if it can support/deny such a correlation. More

In line with the results of present study, we suggest that the exploitation of forage resources by agro-pastoralists occurs in a non-random manner. The use of forage resources is guided by a series of functional characters related to palatability and nutritional value, which determine preferential use due to the better quality of resource. At the same time, we understand that forage uses are complex and multifactorial in nature, and regulated in a substantial way by seasonality and ecological factors (Fig. 5), such as the availability of plant resources and local diversity.Figure 5Diagrammatic representation for the effects of vegetation cover and seasonality on forage resource selection in Dry Forests. Image created with Microsoft Office 2019 PowerPoint (www.office.com).Full size imageThe differences of plant species cited between areas reveal the positive effect of vegetation cover on the use and knowledge of plants by agro-pastoralists. Our findings reveal that the greater number of plant species mentioned by agro-pastoralists in Area II is directly associated with greater availability of resources in this area, as long as we consider vegetation cover as availability of resources, which allows different species to be used throughout the year. On the other hand, in regions with low vegetation cover (Area I), the low availability of resources limits the use and knowledge of plants by residents, which can lead to greater pressure on a small set of available species. Such findings reinforce the importance of vegetation cover for ecosystem provision of goods and services to human populations that depend directly or indirectly on these services.The most represented families found in the present study have also been reported in several other ethnobotanical studies6,16,17,29, with emphasis on Fabaceae and Poaceae, which are recognized for their high forage potential, which derives, above all, from high palatability and nutritional value30. Simultaneously, citations mostly for native species reflect the importance and potential of Caatinga resources as important components of the ruminant diet11, both for the woody and herbaceous strata, corroborating the estimate in the literature that 70% of vegetation has potential use as forage31.The characteristic seasonality of vegetation, on the other hand, represents a limiting factor for forage productivity, culminating in high fluctuations in quality and availability, as well as changes in the dominance of different strata and composition of forage species throughout the seasons11,32. The seasonal distribution of species explains the similarity of seasons between areas, with a higher similarity percentage for the dry seasons, since there is less availability of resources to be exploited compared to the rainy season. In this context, the potentially used species are commonly accessible woody species in both areas. However, during the rainy season, the high availability of herbaceous plants regulates different uses (Fig. 4), but even so, they also exhibit relatively similar patterns, mainly due to the woody component that denotes the common demand by ruminants at the beginning of this season.The effect of climatic variables on vegetation use patterns was documented by16,17, both of which showed greater richness in the use of herbaceous forage during the rainy season, a finding that reflects the seasonal distribution—restriction to that season—and decrease in the qualitative character of annual species33. At the same time, it also reflects the greater number of unique species for the rainy season. However, when compared to woody strata, significant differences in terms of richness are not found because although the diversity of herbaceous species in the Caatinga is greater24, it is much less known than that of the tree-shrub stratum11.Agro-pastoralists even characterize animal preferences for herbaceous stratum, but as its diversity is immense and ephemeral, they claim to have limited ability to identify the species. The high abundance of resources in the rainy season also reduces the concern with forage use, which implies less attention to the species that are consumed. In contrast, woody species, due to multiple uses and greater availability over time, tend to be better known10,34, with a different effect in the dry season making the optimal foraging pattern in this period inherent to the knowledge of agro-pastoralists35.In addition, according to the ecological appearance hypothesis, there is a general tendency for less apparent species to be neglected by populations36. Some studies have corroborated the hypothesis within the context of forage use, with woody species being cited more and having more uses6,15. In addition, people tend to focus on resources whose supply is given continuously10, which may explain why woody species are well represented in both seasons.Security in the provisioning of ecosystem services is an essential component for local populations, and thus woody species are highly valued because they reflect predictability of use15,35. This can be a particularly influential criterion because perennial or late leaf deciduous species, such as Cynophalla flexuosa and Myracrodruon urundeva, had significant amounts of citations and perceptions employing high valuation, as represented by some statements by some interviewees: “É um refrigero na seca” (it is savage in the dry season), “É uma ração boa na seca” (it is a good food in the dry season).In turn, differences in richness of the species cited by the two areas corroborate our first hypothesis that populations inserted in environments with greater vegetation cover tend to cite more species. In line with these findings, considerable floristic dissimilarity was also found between the two areas, given the exclusivity of species. Such dissimilarity may suggest particularities in the vegetation attributes of each area, such as greater floristic diversity7,37,38.Since anthropic processes are irregularly distributed in space, variation in the provisioning of ecosystem services by vegetation also occurs, and influences different collection profiles39. On the other hand, areas with greater species richness have been shown to have greater use patterns6,7. The larger number of species cited as woody and native for Area II is, therefore, associated with greater general richness, as well as herbaceous species present in the rainy season. In contrast, common species are reflected in trends of similar foraging patterns, as well as the presence of common species between areas38. In addition to different levels of disturbance, differences in floristic composition between areas may also be due to edaphic variation40.Our second hypothesis was refuted because the difference in the richness of exotic species between the areas. Plausible explanations for this finding are that, in general, exotic herbaceous species are commonly used for forage in the semi-arid region of Brazil41. Herbaceous species comprise the primary component of the ruminant diet. However, in the midst of their occurrence restricted to the short rainy period, exotic species, mainly of Fabaceae and Poaceae, have been introduced to increase the forage availability, which currently represents an important attribute of forage resources in the Caatinga41,42,43. At the same time, and to also increase the availability of forage resources, the cultivation of species by agro-pastoralists may be common in their properties44, mainly exotics, such as Prosopis juliflora, that have high adaptive potential and governmental incentives45.Regarding use patterns, according to the data presented here it is possible to state that agro-pastoralists ’ experiences with herding activities provide an accumulation of a vast knowledge about forage resources15. This knowledge allows forage resources to be characterized by their potential according to a variety of criteria associated with seasonal variation and qualitative attributes, as commonly found by other studies14,15,16,17,37. Such criteria are often revealed by qualitative approaches that define the valuation perception of resources. Thus, nutritional value and palatability can be implicitly associated with the definitions of “É uma ração boa” (it is a good food), “o bicho gosta muito” (the animals like it very much) and “Rico em proteínas” (rich in protein).It should be added that the establishment of intrinsic relationships with resources allows a particular understanding at a high level of detail15,35, such as changes in palatability throughout development with descriptions including chemical17 and structural changes. Studies confirm that some Caatinga species vary in their chemical composition during leaf maturation, which influences nutritional quality17,46.In addition to revealing the domain of information, this body of knowledge allows maximizing forage use based on nutritional properties weighted by availability14,37. Nunes37 confirmed that the forage species selected by informants and the criteria they adopted coincided with nutritional values measured by the literature, and that, as also found in the present study, younger plants were recognized as highly appreciated by animals. This appreciation is due to the greater palatability of plant organs at this stage47. This is a matter of concern for the sustainability of the Caatinga, since direct or indirect grazing has compromised the regeneration process12 since younger individuals are clearly more sensitive to damage48.Also, considering the potential of Caatinga, we suggest that investment through government actions encourage the cultivation of native species to ensure the production of forage and, consequently, guarantee the sustainability of livestock activity and the ecosystem in question. More