Ecology

At least 11 coastal ecosystems have been considered based on a minimum of actionably defined criteria to be “blue carbon ecosystems”. These include mangrove wetlands, tidal marshes (salt, brackish, fresh), seagrasses, salt flats, freshwater (upper estuarine) tidal forests, macroalgae, phytoplankton, coral reef, marine fauna (fish), oyster reefs, and mud flats5; all but three of these would be considered wetlands, with salt flats and mud flats being examples of non-emergent (plant) blue carbon wetlands. Herein, we focus on mangroves.Adjusting intrinsic leaf-level photosynthetic water use efficiency (({WUE}_{int})) in response to environmental gradients (Introduction)We used data provided by B.F. Clough & R.G. Sims20, which presented leaf-scale net photosynthesis (({P}_{n}) [sic]; μmol CO2 m−2 s−1), stomatal conductance (({g}_{w}): mol m−2 s−1), leaf-intercellular CO2 concentrations (({c}_{i}): μl l−1), and intrinsic photosynthetic water use efficiency (({WUE}_{int}): (frac{{P}_{n}}{{g}_{w}})) for 19 mangrove species occupying 9 different sites in Papua New Guinea and northern Australia. These field data were collected using an infrared gas analyzer (model Li-6000, Li-Cor Biosciences, Inc., Lincoln, NE, USA) attached to leaves at saturation light levels (reported as  > 800 μmol PPFD m−2 s−1). Soil salinity at the time of data collection ranged from 10 to 49 psu, and median long-term atmospheric temperature and relative humidity among sites ranged from 19.9 to 27.4 °C and 35.1 to 92.2%, respectively (Fig. S1)20. These data were among the first to offer insight from field study into the plasticity of mangroves across a range of natural salinity and aridity gradients to adjust leaf-level ({WUE}_{int}) as needed for local environmental condition. While it is not new for trees to adjust their ({WUE}_{int}) when they develop in arid, semi-arid, or even some humid and tropical environments60, what is distinctive is that mangroves may be further driven to water savings by salinity gradients as a condition of development.
({{varvec{W}}{varvec{U}}{varvec{E}}}_{{varvec{i}}{varvec{n}}{varvec{t}}}) and individual tree water use of mangrove wetlands versus terrestrial ecosystemsFor Fig. 1a, we compare leaf-level ({WUE}_{int}) data collected from 17 published papers (using maximum and minimum values), providing 67 independent measurements of ({WUE}_{int}) for mangroves (Table S3). While we mention in the main text that as many as 214 independent measurements of water use efficiency are available, not all of these present raw ({P}_{n}) or ({g}_{w}) data, with some reporting leaf transpiration (({T}_{r})) which do not enable reporting of intrinsic water use efficiencies. Also, we strategically included studies from reproducible experimental designs and readily available papers. Mangrove species included in this review represented a global distribution of greenhouse and field observations, and encompassed species in the following mangrove genera: Rhizophora, Avicennia, Laguncularia, Bruguiera, Aegialitis, Aegiceras, Ceriops, Sonneratia, Kandelia, Excoecaria, Heritiera, Xylocarpus, and Conocarpus.We then accessed an existing database (n = 11,328 observations) that reported raw ({P}_{n}) and ({g}_{w}) data from 210 upland deciduous and evergreen shrubs and trees of savannah, boreal, temperate, and tropical habitats60. From these data, we evaluated a range of upland tree and shrub species that occurred and developed naturally in environments along a global gradient of vapor pressure deficit (i.e., atmospheric moisture and temperature), including arid, semi-arid, dry semi-humid, and humid locations.For Fig. 1b, we started with a review by Wullschleger et al.61 that provides maximum individual tree water use data (L H2O day−1) from 52 published studies representing 67 species of upland trees from around the world. Of those studies, dbh values (8 to 134 cm) were provided alongside 47 individual tree water use values. Maximum individual tree water use and dbh (4.1 to 45.3 cm) were available from the original source for 8 mangrove studies representing 7 species from French Guiana, Mayotte Island (Indian Ocean), China, Florida (USA), and Louisiana (USA) (Table S4). These represent the extent of published sap flow data that provided both individual tree water use and dbh from mangroves (numerically); e.g., we could not extract specific individual tree water use versus dbh from a Moreton Bay (Australia) study site62, south Florida study site63, or from five additional study sites in China51,52. However, regressions for two of the Chinese study sites provided over two years51 indicated that mangrove trees from a suite of species ranging in dbh from 8 to 24 cm used approximately 0.76 and 9.31 L H2O day−1, or 0.53 L H2O day−1 cm−1 of dbh. These apparent rates were even lower than what was reported as average for mangroves in Fig. 1b of 1.4 L H2O day−1 cm−1. The mangrove species included in this analysis were Avicennia germinans (L.) L., Laguncularia racemosa (L.) C.F. Gaertn., Rhizophora mangle L., Ceriops tagal (Perr.) C.B. Rob., Rhizophora mucronata Lam., Sonneratia apetala Buch.-Ham, and Sonneratia caseolaris (L.) Engl.. Additional comparative mangrove species reported by B. Leng & K.-F. Cao51 included Bruguiera sexangula (Lour.) Poir., Bruguiera sexangula var. rhynchopetala W.C. Ko, Excoecaria agallocha L., Rhizophora apiculata Blume, Sonneratia alba Sm., and Xylocarpus granatum J. Koenig.Estimation of canopy transpiration (({{varvec{E}}}_{{varvec{c}}})) from net primary productivity dataEstimates of carbon uptake from CO2 can provide insight into the water use requirement for that uptake of carbon64. We used leaf-level instantaneous water use efficiency (({WUE}_{ins}): (frac{{P}_{N}}{{T}_{r}})), which relates to net CO2 uptake from leaves of the dominant mangrove forest canopy relative to the specific amount of water used, and developed a predictive relationship (predicted) for determining mangrove net primary productivity (NPP) values from ({E}_{c}) using ({WUE}_{ins}). For A. germinans, L. racemosa, and R. mangle forest components, we used light-saturated, leaf-level ({WUE}_{ins}) values of 3.82 ± 0.3, 4.57 ± 0.3, and 5.15 ± 0.4 mmol CO2 (mol H2O)−1 [± 1 SE], respectively, from mangrove saplings and small trees of south Florida65. ({WUE}_{ins}) values were stratified by species relative to basal area distributions on each south Florida study plot, converted from molar fractions of H20 (from ({E}_{c}) determination) and CO2 to molecular weights, and multiplied by ({WUE}_{ins}) with applicable unit conversions to attain kg CO2 m−2 year−1. This value was multiplied by 0.273 to yield kg C m−2 year−1.This predictive relationship was validated in two independent ways. First, for one of the calibration sites (lower Shark River, Everglades National Park, Florida, USA), we modeled ({E}_{c}) from sap flow data50, determined NPP from ({WUE}_{ins}) calculations relative to the amount of water the stand used, and had independent measurements of net ecosystem exchange (NEE) of CO2 between the mangrove ecosystem and atmosphere from an eddy flux tower66. For this site, NPP estimation and NEE were closely aligned once soil CO2 effluxes were accounted; respiratory CO2 effluxes from soil and pneumatophores were determined to be 1.2 kg C m−2 year−1 from previous study67. Using our NPP estimations from ({WUE}_{ins}) calculations and subtracting soil and pneumatophore CO2 effluxes of 1.2 kg C m−2 for 2004 and 0.8 kg C m−2 for 2005 (partial year), NPP becomes 0.96 kg C m−2 for 2004 and 0.85 kg C m−2 for January to August of 2005 (see Observed 1, Florida on Fig. S2). Our approach underestimated NPP from ({E}_{c}) relative to measurements from eddy covariance by 0.21 kg C m-2 for 2004 (within 17.5% of predicted) and was nearly identical for 2005 (within 0.02 kg C m−2, or 2% of predicted).Second, we wanted to determine whether ({E}_{c})-to-NPP predictions developed on a few sites in south Florida, USA, represented other global sites, so we included an analysis from several mangrove sites in Guangdong Province, China, to represent an entirely different location. Similar to south Florida analyses, we combined data for NPP from co-located sites of ({E}_{c}) determination using sap flow techniques. NPP of the mangrove forests were measured using multiple procedures (including eddy flux) for improved accuracy68,69. The relationships of predicted NPP versus ({E}_{c}) and observed NPP versus ({E}_{c}) did not differ between Florida and China (t = 0.48, p = 0.643).Projecting mangrove ({{varvec{E}}}_{{varvec{c}}}) to other locationsWe reviewed data from 26 published records that report mangrove NPP, or enough data to estimate NPP, from 71 study sites located in the Florida-Caribbean Region (25 sites) and Asia–Pacific Region (46 sites) (Table S5). Table S1 reveals mangrove literature sources used, as well as how NPP was estimated from values provided in the original source; itemizes assumptions for determinations of aboveground NPP, wood production, litter production, and root production from various ratios70; and reveals unit conversions.We then convert NPP to ({E}_{c}) for all 71 sites using the predicted curve in Fig. S2 (Eq. 1, main text), and provide summary results by location in Table S1. Regional (ET) data were extracted from the MODIS Global Evapotranspiration Project (MOD16-A3), which are provided at a resolution of 1-km. The locations of mangrove NPP study sites were identified, assigned to a single 1-km2 grid in MOD16, and (ET) was extracted from that grid and used for ({E}_{c})-to-(ET) comparison. Average (ET) from single cells (1 km2) was combined with the average of up to 8 additional neighboring cells to provide comparative (ET) projections over up to 9 km2 for each location from 2000 to 2013 to compare sensitivity among suites of the specific MODIS16-A3 cells selected over land. When neighboring cells were completely over water, they were excluded since component mangrove forest ({E}_{c}) estimation was not possible from the cells. Estimates of (ET) by individual cells used to compare with mangrove ({E}_{c}) versus up to 9 cells differed by an average of only 43 mm H2O year−1 (± 16 mm H2O year−1, S.E.). Therefore, we use (ET) from individual, overlapping ({E}_{c}) cells in Table S1.The average ({E}_{c})-to-(ET) ratio from mangroves was subtracted from ({E}_{c})-to-(ET) ratio for specific ecoregions48, and this ratio difference was assumed to represent net water use strategy affecting differences by the dominant vegetation between ecosystem types. We were also mindful that salinity reductions can affect ({E}_{c}). We used scaled (0–1) mean and standard deviations from ({WUE}_{int}) data previously reported for mangroves (Fig. S1)20. Standard deviation was 32% of mean ({WUE}_{int}) related to salinity gradients, and if we re-scale this deviation to ({E}_{c}) data and add it to the mean ({E}_{c}) to assume low salinity, average ({E}_{c})-to-(ET) ratio becomes 57.4%. This is theoretical and assumes a relatively linear relationship between ({WUE}_{int}) and ({E}_{c}).Comparative water use scaling among ecoregionsTable 1 presents the projected reduction in water used through ({E}_{c}) if a mangrove ({E}_{c})-to-(ET) ratio was applied to tropical rainforest (290.52 mm H2O year−1), temperate deciduous forest (131.76 mm H2O year−1), tropical grassland (110.77 mm H2O year−1), temperate grassland (46.48 mm H2O year−1), temperate coniferous forest (54.96 mm H2O year−1), desert (22.99 mm H2O year−1), and Mediterranean shrubland (12.08 mm H2O year−1). To convert potential water use differences to kL H2O ha−1 year−1 (as presented in the abstract), the following calculation is used (using the example of tropical rainforest):$$frac{290.52 L {H}_{2}O {year}^{-1}}{1 {m}^{2}} times frac{mathrm{10,000 }{m}^{2}}{1 ha} times frac{1 kL {H}_{2}O}{mathrm{1000 }L {H}_{2}O} =mathrm{2905 } kL {H}_{2}O {ha}^{-1}{year}^{-1}$$
(2)
For comparisons made to mature ( > 12 years) oil palm (Elaeis guineensis Jacq.) plantations, ({E}_{c})-to-(ET) ratio was assumed to range from 5332 to 70%33, for a water use difference of 1170 and 3160 kL H2O ha−1 year−1, respectively, relative to annual global mangrove (ET) (of 1172 mm). We multiply these values by the 18,467 ha of land area that was converted from mangroves to oil palm31 to attain potential water use differences of 21.6 to 58.4 GL H2O year−1 from avoided conversion of mangrove to oil palm in this region.Global water use scalingIn order to determine how much global mangrove area is adjacent to each ecoregion, we conducted a cross-walk between terrestrial ecoregions71 and those used by Global Mangrove Watch in the 2010 classification of global mangrove area72. Terrestrial ecoregions used by Schlesinger & Jasechko48 were then able to be associated with specific mangrove areas (Table S6). In other words, given a specific ecoregion, we determined how much mangrove area would be occurring within that same ecoregional geography. Global mangrove area assignment to those ecoregions mapped within 0.1% of the total mangrove area of 13,760,000 ha reported in Bunting et al.72. To convert kL H2O ha−1 year−1 to GL H2O year−1 among ecoregions, the following calculation was used (continuing with the example of tropical rainforest, which has an area of adjacent mangroves of 112,331.9 km2):$$frac{mathrm{2905} kL {H}_{2}O {year}^{-1}}{1 ha} times frac{100 ha}{1 {km}^{2}} times frac{mathrm{112,331.90} {km}^{2 }mangroves}{1.0 times {10}^{6} kL {H}_{2}0} times frac{1 GL {H}_{2}O}{1} = mathrm{32,632.42} GL {H}_{2}O {year}^{-1}$$
(3)
Agent-based modelling of individual tree water use (Discussion)The BETTINA model simulates the growth of mangrove trees as a response to above- and below-ground resources, i.e. light and water41. In the model, an individual tree is described by four geometric measures, including stem radius, stem height, crown radius and root radius; attributing functional relevance in terms of resource uptake. Aiming to maximize resource uptake, new biomass is allocated to increase these measures in an optimal but not constant proportion. Water uptake of the tree is driven by the water potential gradient between the soil and the leaves. Thus, porewater salinity is part of what determines the water availability for plants.With the BETTINA model, we simulated the growth of nine individual mangrove trees under different salinity conditions, ranging between 0 and 80 psu, while all other environmental and tree-specific conditions were kept constant. Simulation time was 200 years so that trees could achieve very close to their maximum possible size, and the hydrological parameters were similar to that reported previously42. We can show that the ratio of the actual transpiration to the potential transpiration decreases with increasing salinity; plants use less water. Potential transpiration was the transpiration of a given tree without a simulated reduction in water availability due to porewater salinity. These parameter details are presented graphically for mangroves (Fig. S3), comparing porewater salinity along a gradient against the ratio of actual-to-potential individual tree transpiration.Further, BETTINA simulation results include morphological plasticity adjustments to allometry. To highlight this, we also displayed results assuming a constant allometry as for 40 psu. Naturally, for this arbitrary benchmark the solid and the dashed line coincide (Fig. 3a). Adaptation to higher salinities improves water uptake (primarily girth and root growth), thus the adapted trees (solid lines) have a higher water uptake than the average allometry (dashed lines) for salinities below 40 psu. Lower salinities promote increase of height and crown radius to improve light availability. That is why the adapted trees have a lower water uptake than an average tree would for salinities above 40 psu. Tree water use decreases with increasing salinity (Fig. 3b), as ({WUE}_{int}) coincidently increases (Fig. S1).Virtual water use explained (Discussion)Water is required to produce products or acquire services from natural ecosystems; e.g., forest products, fisheries biomass, nutrient processing (nitrification, denitrification), food production. If a net kilogram of a food is grown on a hectare of land where water is abundant and that kilogram of food requires 400 mm of water to be produced, the export of that food to an area of low water availability provides an ecosystem service in the amount of 1 kg of food, plus 400 mm of “virtual” water not actually needed at the destination but used at the source. This water is defined as the product’s “virtual water content”56. There is a rich body of literature exploring the concept of virtual water73,74, but we expand on this concept here as a comparison among 7 ecoregions48 and mangroves. Raw data used for calculations are presented in Table S2.Statistical analysisData for leaf-level ({WUE}_{int}) comparisons between terrestrial woody plants and mangroves, as well as individual tree water use by dbh for both terrestrial and mangrove trees, were not normally distributed. We used a Kruskal–Wallis ANOVA based on ranks, and the Dunn’s Method for difference tests. Individual tree water use by dbh for both terrestrial and mangrove trees were determined using linear regression, mostly applied to mean values. For a couple of mangrove studies, only median values were extractable from minimum and maximum values. Likewise, all other data relationships were best fit with linear models, including the calibration curves between ({E}_{c}) and NPP. All data were analyzed using SigmaPlot (v. 14.0, Systat, Inc., Palo Alto, California, USA). More

Hoham, R. W. & Remias, D. Snow and glacial algae: A review. J. Phycol. 56, 264–282 (2020).PubMed 
PubMed Central 

Google Scholar 
Williamson, C. J. et al. Glacier algae: A dark past and a darker future. Front. Microbiol. 10, 524 (2019).
PubMed 
PubMed Central 

Google Scholar 
Hodson, A. J., Mumford, P. N., Kohler, J. & Wynn, P. M. The High Arctic glacial ecosystem: New insights from nutrient budgets. Biogeochemistry 72, 233–256 (2005).CAS 

Google Scholar 
Stibal, M., Tranter, M., Telling, J. & Benning, L. G. Speciation, phase association and potential bioavailability of phosphorus on a Svalbard glacier. Biogeochemistry 90, 1–13 (2008).CAS 

Google Scholar 
Telling, J. et al. Microbial nitrogen cycling on the Greenland ice sheet. Biogeosciences 9, 2431–2442 (2012).ADS 
CAS 

Google Scholar 
Cook, J. M. et al. Glacier algae accelerate melt rates on the south-western Greenland ice sheet. Cryosphere 14, 309–330 (2020).ADS 

Google Scholar 
Yallop, M. L. et al. Photophysiology and albedo-changing potential of the ice algal community on the surface of the Greenland ice sheet. ISME J. 6, 2302–2313 (2012).CAS 
PubMed 
PubMed Central 

Google Scholar 
Stibal, M. et al. Algae drive enhanced darkening of bare ice on the Greenland ice sheet. Geophys. Res. Lett. 44, 11463–11471 (2017).ADS 

Google Scholar 
Di Mauro, B. et al. Glacier algae foster ice-albedo feedback in the European Alps. Sci. Rep. 10, 1–9 (2020).ADS 

Google Scholar 
Lutz, S. et al. The biogeography of red snow microbiomes and their role in melting arctic glaciers. Nat. Commun. 7, 1–9 (2016).ADS 

Google Scholar 
Ganey, G. Q., Loso, M. G., Burgess, A. B. & Dial, R. J. The role of microbes in snowmelt and radiative forcing on an Alaskan icefield. Nat. Geosci. 10, 754–759 (2017).ADS 
CAS 

Google Scholar 
Khan, A. L., Dierssen, H. M., Scambos, T. A., Höfer, J. & Cordero, R. R. Spectral characterization, radiative forcing and pigment content of coastal Antarctic snow algae: Approaches to spectrally discriminate red and green communities and their impact on snowmelt. Cryosphere 15(1), 133–148 (2021).ADS 

Google Scholar 
Huovinen, P., Ramírez, J. & Gómez, I. Remote sensing of albedo-reducing snow algae and impurities in the Maritime Antarctica. ISPRS J. Photogramm. Remote Sens. 146, 507–517 (2018).ADS 

Google Scholar 
Stibal, M., Šabacká, M. & Žárský, J. Biological processes on glacier and ice sheet surfaces. Nat. Geosci. 5, 771–774 (2012).ADS 
CAS 

Google Scholar 
Williamson, C. J. et al. Algal photophysiology drives darkening and melt of the Greenland ice sheet. Proc. Natl. Acad. Sci. https://doi.org/10.1073/pnas.1918412117 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
Chevrollier, L.-A. et al. Light absorption and albedo reduction by pigmented microalgae on snow and ice. J. Glaciol. https://doi.org/10.1017/jog.2022.64 (2022).Article 

Google Scholar 
Procházková, L., Řezanka, T., Nedbalová, L. & Remias, D. Unicellular versus filamentous: The glacial alga ancylonema alaskana comb. et. stat. nov. and its ecophysiological relatedness to ancylonema nordenskioeldii (zygnematophyceae, streptophyta). Microorganisms 9(5), 1103 (2021).PubMed 
PubMed Central 

Google Scholar 
Remias, D. et al. Characterization of an UV- and VIS-absorbing, purpurogallin-derived secondary pigment new to algae and highly abundant in Mesotaenium berggrenii (Zygnematophyceae, Chlorophyta), an extremophyte living on glaciers. FEMS Microbiol. Ecol. 79, 638–648 (2012).CAS 
PubMed 

Google Scholar 
Williamson, C. J. et al. Ice algal bloom development on the surface of the Greenland ice sheet. FEMS Microbiol. Ecol. 94, 1–10 (2018).
Google Scholar 
Remias, D., Lütz-Meindl, U. & Lütz, C. Photosynthesis, pigments and ultrastructure of the alpine snow alga Chlamydomonas nivalis. Eur. J. Phycol. 40, 259–268 (2005).CAS 

Google Scholar 
Leya, T., Rahn, A., Lütz, C. & Remias, D. Response of arctic snow and permafrost algae to high light and nitrogen stress by changes in pigment composition and applied aspects for biotechnology. FEMS Microbiol. Ecol. 67, 432–443 (2009).CAS 
PubMed 

Google Scholar 
Müller, T., Bleiß, W., Martin, C. D., Rogaschewski, S. & Fuhr, G. Snow algae from northwest Svalbard: Their identification, distribution, pigment and nutrient content. Polar Biol. 20, 14–32 (1998).
Google Scholar 
Remias, D., Karsten, U., Lütz, C. & Leya, T. Physiological and morphological processes in the Alpine snow alga Chloromonas nivalis (Chlorophyceae) during cyst formation. Protoplasma 243, 73–86 (2010).PubMed 

Google Scholar 
Bidigare, R. R. et al. Evidence for a photoprotective function for secondary carotenoids of snow algae taxonomy, life histories, ecology and geographical habitats and polar regions. J. Phycol. 434, 427–434 (1993).CAS 

Google Scholar 
Remias, D. & Lütz, C. Characterisation of esterified secondary carotenoids and of their isomers in green algae: A HPLC approach. Arch. Hydrobiol. Suppl. Algol. Stud. 124, 85–94 (2007).CAS 

Google Scholar 
Dial, R. J., Ganey, G. Q. & Skiles, S. M. What color should glacier algae be ? An ecological role for red carbon in the cryosphere. FEMS Microbiol. Ecol. 94(3), 1–9. https://doi.org/10.1093/femsec/fiy007 (2018).Article 
CAS 

Google Scholar 
Uusikivi, J., Vähätalo, A. V., Granskog, M. A. & Sommaruga, R. Contribution of mycosporine-like amino acids and colored dissolved and particulate matter to sea ice optical properties and ultraviolet attenuation. Limnol. Oceanogr. 55, 703–713 (2010).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Xie, H., Aubry, C., Zhang, Y. & Song, G. Chromophoric dissolved organic matter (CDOM) in first-year sea ice in the western Canadian Arctic. Mar. Chem. 165, 25–35 (2014).CAS 

Google Scholar 
Holzinger, A. & Lütz, C. Algae and UV irradiation: Effects on ultrastructure and related metabolic functions. Micron 37, 190–207 (2006).PubMed 

Google Scholar 
Piiparinen, J. et al. The contribution of mycosporine-like amino acids, chromophoric dissolved organic matter and particles to the UV protection of sea-ice organisms in the Baltic Sea. Photochem. Photobiol. Sci. 14, 1025–1038 (2015).CAS 
PubMed 

Google Scholar 
Cook, J. M. et al. Quantifying bioalbedo: A new physically based model and discussion of empirical methods for characterising biological influence on ice and snow albedo. Cryosphere 11, 2611–2632 (2017).ADS 

Google Scholar 
Lutz, S., Anesio, A. M., Jorge Villar, S. E. & Benning, L. G. Variations of algal communities cause darkening of a Greenland glacier. FEMS Microbiol. Ecol. 89, 402–414 (2014).CAS 
PubMed 

Google Scholar 
Hoham, R. W. & Ling, H. U. Snow algae: Tile effects of chemical and physical factors on their life cycles and populations. In Journey to Diverse Microbial Worlds (ed. Seckback, J.) 131–145 (Springer, Netherlands, 2000).
Google Scholar 
Procházková, L., Leya, T., Krížková, H. & Nedbalová, L. Sanguina nivaloides and Sanguina aurantia gen. Et spp. Nov. (Chlorophyta): The taxonomy, phylogeny, biogeography and ecology of two newly recognised algae causing red and orange snow. FEMS Microbiol. Ecol. 95, 1–21 (2019).
Google Scholar 
Hillebrand, H., Dürselen, C. D., Kirschtel, D., Pollingher, U. & Zohary, T. Biovolume calculation for pelagic and benthic microalgae. J. Phycol. 35, 403–424 (1999).
Google Scholar 
Montagnes, D. J. S., Berges, J. A., Harrison, P. J. & Taylor, F. J. R. Estimating carbon, nitrogen, protein and chlorophyll a from volume in marine phytoplankton. Limnol. Oceanogr. 39, 1044–1060 (1994).ADS 
CAS 

Google Scholar 
Schreiber, U., Gademann, R., Ralph, P. J. & Larkum, A. W. D. Assessment of photosynthetic performance of Prochloron in Lissoclinum patella in hospite by chlorophyll fluorescence measurements. Plant Cell Physiol. 38, 945–951 (1997).CAS 

Google Scholar 
Jassby, A. D. & Platt, T. Mathematical formulation of the relationship between photosynthesis and light for phytoplankton. Limnol. Oceanogr. 21, 540–547 (1976).ADS 
CAS 

Google Scholar 
Silsbe, G. M. & Malkin, S. Y. Package ‘phytotools’: Phytoplankton Production Tools. (2015).Aigner, S., Remias, D., Karsten, U. & Holzinger, A. Unusual phenolic compounds contribute to ecophysiological performance in the purple-colored green alga Zygogonium ericetorum (Zygnematophyceae, Streptophyta) from a high-alpine habitat. J. Phycol. 49, 648–660 (2013).CAS 
PubMed 
PubMed Central 

Google Scholar 
Holzinger, A. et al. Arctic, Antarctic, and temperate green algae Zygnema spp. under UV-B stress: Vegetative cells perform better than pre-akinetes. Protoplasma 255, 1239–1252 (2018).CAS 
PubMed 
PubMed Central 

Google Scholar 
Bidigare, R. R., Ondrusek, M. E., Morrow, J. H. & Kiefer, D. A. In-vivo absorption properties of algal pigments. Ocean Opt. X 1302, 290 (1990).ADS 

Google Scholar 
Clementson, L. A. & Wojtasiewicz, B. Dataset on the absorption characteristics of extracted phytoplankton pigments. Data Br. 24, 103875 (2019).
Google Scholar 
Bidigare, R. R., Ondrusek, M. E., Morrow, J. H. & Kiefer, D. A. In vivo absorption properties of algal pigments. Ocean Opt. X 1302, 290–302 (1990).ADS 

Google Scholar 
Bricaud, A. & Stramski, D. Spectral absorption coefficients of living phytoplankton and nonalgal biogenous matter: A comparison between the Peru upwelling areaand the Sargasso Sea. Limnol. Oceanogr. 35, 562–582 (1990).ADS 
CAS 

Google Scholar 
Duval, B., Shetty, K. & Thomas, W. H. Phenolic compounds and antioxidant properties in the snow alga Chlamydomonas nivalis after exposure to UV light. J. Appl. Phycol. 11, 559–566 (2000).
Google Scholar 
Onuma, Y. et al. Observations and modelling of algal growth on a snowpack in north-western Greenland. Cryosphere 12, 2147–2158 (2018).ADS 

Google Scholar 
Christner, B. C. et al. Microbial processes in the weathering crust aquifer of a temperate glacier. Cryosphere 12, 3653–3669 (2018).ADS 

Google Scholar 
Cook, J. M., Hodson, A. J. & Irvine-Fynn, T. D. L. Supraglacial weathering crust dynamics inferred from cryoconite hole hydrology. Hydrol. Process. 30, 433–446 (2016).ADS 

Google Scholar 
Tedstone, A. J. et al. Algal growth and weathering crust state drive variability in western Greenland ice sheet ice albedo. Cryosphere 14, 521–538 (2020).ADS 

Google Scholar 
Smith, L. C. et al. Efficient meltwater drainage through supraglacial streams and rivers on the southwest Greenland ice sheet. Proc. Natl. Acad. Sci. U. S. A. 112, 1001–1006 (2015).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Irvine-Fynn, T. D. L. et al. Storage and export of microbial biomass across the western Greenland ice sheet. Nat. Commun. 12, 1–11 (2021).
Google Scholar 
Cameron, K. A. et al. Meltwater export of prokaryotic cells from the Greenland ice sheet. Environ. Microbiol. 19, 524–534 (2017).PubMed 

Google Scholar 
Stibal, M. et al. Environmental controls on microbial abundance and activity on the Greenland ice sheet: A multivariate analysis approach. Microb. Ecol. 63, 74–84 (2012).PubMed 

Google Scholar 
Mernild, S. H., Liston, G. E., Hasholt, B. & Knudsen, N. T. Snow distribution and melt modeling for Mittivakkat Glacier, Ammassalik Island, southeast Greenland. J. Hydrometeorol. 7, 808–824 (2006).ADS 

Google Scholar 
Stibal, M., Elster, J., Šabacká, M. & Kaštovská, K. Seasonal and diel changes in photosynthetic activity of the snow alga Chlamydomonas nivalis (Chlorophyceae) from Svalbard determined by pulse amplitude modulation fluorometry. FEMS Microbiol. Ecol. 59, 265–273 (2007).CAS 
PubMed 

Google Scholar 
Remias, D., Holzinger, A., Aigner, S. & Lu, C. Ecophysiology and ultrastructure of Ancylonema nordenskioldii (Zygnematales, Streptophyta), causing brown ice on glaciers in Svalbard (high arctic). 899–908 (2012). doi:https://doi.org/10.1007/s00300-011-1135-6Remias, D., Holzinger, A. & Lütz, C. Physiology, ultrastructure and habitat of the ice Alga Mesotaenium berggrenii (Zygnemaphyceae, Chlorophyta) from Glaciers in the European Alps. Phycologia 48, 302–312 (2009).
Google Scholar 
Nakashima, T. et al. Spatial and temporal variations in pigment and species compositions of snow algae on Mt. Tateyama in Toyama prefecture, Japan.. Front. Plant Sci. 12, 1–16 (2021).
Google Scholar 
Remias, D., Albert, A. & Lütz, C. Effects of realistically simulated, elevated UV irradiation on photosynthesis and pigment composition of the alpine snow alga Chlamydomonas nivalis and the arctic soil alga Tetracystis sp. (Chlorophyceae). Photosynthetica 48, 269–277 (2010).CAS 

Google Scholar 
Procházková, L., Remias, D., Holzinger, A., Řezanka, T. & Nedbalová, L. Ecophysiological and ultrastructural characterisation of the circumpolar orange snow alga Sanguina aurantia compared to the cosmopolitan red snow alga Sanguina nivaloides (Chlorophyta). Polar Biol. 44, 105–117 (2021).PubMed 

Google Scholar 
Niyogi, K. K. Photoprotection revisited: Genetic and molecular approaches. Annu. Rev. Plant Biol. 50, 333–359 (1999).CAS 

Google Scholar 
Sakshaug, E. & Holm-hansen, O. Photoadaptation in Antarctic phytopfankton: Variations in growth rate, chemical composition and P versus I curves. J. Plankton Res. 8, 459–473 (1986).
Google Scholar 
Malerba, M. E., Palacios, M. M., Palacios Delgado, Y. M., Beardall, J. & Marshall, D. J. Cell size, photosynthesis and the package effect: An artificial selection approach. New Phytol. 219, 449–461 (2018).CAS 
PubMed 

Google Scholar 
Wagner, B., Ochs, D. & Bieler, K. Derivatives as antimicrobial agents. Engineering 8, 240–244 (2011).
Google Scholar 
Perini, L. et al. Darkening of the Greenland ice sheet: Fungal abundance and diversity are associated with algal bloom. Front. Microbiol. 10, 557 (2019).PubMed 
PubMed Central 

Google Scholar 
Perini, L. et al. Interactions of fungi and algae from the Greenland ice sheet. Microb. Ecol. https://doi.org/10.1007/s00248-022-02033-5 (2022).Article 
PubMed 

Google Scholar 
Taskjelle, T. et al. Spectral albedo and transmittance of thin young Arctic sea ice. J. Geophys. Res. Ocean. 121, 540–553 (2015).ADS 

Google Scholar 
Lutz, S., Anesio, A. M., Edwards, A. & Benning, L. G. Linking microbial diversity and functionality of arctic glacial surface habitats. Environ. Microbiol. 19, 551–565 (2017).CAS 
PubMed 

Google Scholar 
Smith, H., Dieser, M., McKnight, D., SanClements, M. & Foreman, C. Relationship between dissolved organic matter quality and microbial community composition across polar glacial environments. FEMS Microbiol. Ecol. 94(7), fiy090 (2018).CAS 

Google Scholar 
Kirk, J. T. O. A theoretical analysis of the contributino of algal cells to the attenuation of light within natural waters. New Phytol. 77, 341–358 (1976).
Google Scholar 
Morel, A. & Bricaud, A. Theoretical results concerning light absorption in a discrete medium, and application to specific absorption of phytoplankton. Deep Sea Res. Part A. Oceanogr. Res. Pap. 28, 1375–1393 (1981).ADS 

Google Scholar 
Stuart, V., Sathyendranath, S., Platt, T., Maass, H. & Irwin, B. D. Pigments and species composition of natural phytoplankton populations: Effect on the absorption spectra. J. Plankton Res. 20, 187–217 (1998).CAS 

Google Scholar 
Kirk, J. Light and Photosynthesis in Aquatic Environment (University Press, 1983).
Google Scholar 
Nelson, N. B., Prezelin, B. B. & Bidigare, R. R. Phytoplankton light absorption and the package effect in California coastal waters. Mar. Ecol. Prog. Ser. 94, 217–227 (1993).ADS 

Google Scholar 
Holzinger, A., Allen, M. C. & Deheyn, D. D. Hyperspectral imaging of snow algae and green algae from aeroterrestrial habitats. J. Photochem. Photobiol. B Biol. https://doi.org/10.1016/j.jphotobiol.2016.07.001 (2016).Article 

Google Scholar 
Anesio, A. M. et al. Monitoring glacial algae and impurities on the Greenland Ice Sheet. Aarhus Univ. DCE – Danish Cent. Environ. Energy, Sci. Rep. No. 489 26 (2022).QGIS.org, %Y. QGIS Geographic Information System. QGIS Association. http://www.qgis.org. More

Gibson, R. N. & Yoshiyama, R. M. in Intertidal fishes; life in two worlds (eds. Horn, M. H., Martin, K. L. M. & Chotkowski, M. A.) Intertidal fish communities (Academic Press, 1999).Cox, T. E., Baumgartner, E., Philippoff, J. & Boyle, K. S. Spatial and vertical patterns in the tidepool fish assemblage on the island of O’ahu. Environ. Biol. Fishes 90, 329–342 (2011).Article 

Google Scholar 
Bezerra, L. A. V., Padial, A. A., Mariano, F. B., Garcez, D. S. & Sánchez-Botero, J. I. Fish diversity in tidepools: assembling effects of environmental heterogeneity. Environ. Biol. Fishes 100, 551–563 (2017).Article 

Google Scholar 
Metaxas, A. & Scheibling, R. E. Community structure and organization of tidepools. Mar. Ecol. Prog. Ser. 98, 187–198 (1993).Article 
ADS 

Google Scholar 
Castellanos-Galindo, G. A., Giraldo, A. & Rubio, E. A. Community structure of an assemblage of tidepool fishes on a tropical eastern Pacific rocky shore, Colombia. J. Fish. Biol. 67, 392–408 (2005).Article 

Google Scholar 
Castellanos-Galindo, G. A. & Giraldo, A. Food resource use in a tropical eastern Pacific tidepool fish assemblage. Mar. Biol. 153, 1023–1035 (2008).Article 

Google Scholar 
Moring, J. R. Seasonal absence of fishes in tidepools of a boreal environment (Maine, USA). Hydrobiologia 194, 163–168 (1990).Article 

Google Scholar 
Davis, J. L. D. Spatial and seasonal patterns of habitat partitioning in a guild of southern California tidepool fishes. Mar. Ecol. Prog. Ser. 196, 253–268 (2000).Article 
ADS 

Google Scholar 
Arakaki, S., Tsuchiya, M. & Tokeshi, M. Testing latitudinal patterns of tidepool fish assemblages: local substrate characteristics affect regional-scale trends. Hydrobiologia 733, 45–62 (2014).Article 

Google Scholar 
Arakaki, S. & Tokeshi, M. Assessment of decadal changes in the tidepool fish assemblage of Danjo Islands in the northern East China Sea. Mar. Biodivers. 52, 25 (2022).Article 

Google Scholar 
Bonaca, M. O. & Lipej, L. Factors affecting habitat occupancy of fish assemblage in the Gulf of Trieste (Northern Adriatic Sea). Mar. Ecol. 6, 42–53 (2005).Article 
ADS 

Google Scholar 
Ribeiro, J., Carvalho, G. M., Gonçalves, J. M. S. & Erzini, K. Fish assemblages of shallow intertidal habitats of the Ria Formosa lagoon (South Portugal): influence of habitat and season. Mar Ecol Prog Ser 446, 259–273 (2012).Article 
ADS 

Google Scholar 
Compaire, J. C. et al. Micro- and macroscale factors affecting fish assemblage structure in the rocky intertidal zone. Mar Ecol Prog Ser 610, 175–189 (2019).Article 
ADS 

Google Scholar 
Griffiths, S. P. Spatial and temporal dynamics of temperate Australian rockpool ichthyofaunas. Mar. Freshwater Res. 54, 163–176 (2003).Article 
ADS 

Google Scholar 
Henriques, S., Pais, M. P., Costa, M. J. & Cabral, H. N. Seasonal variability of rocky reef fish assemblages: Detecting functional and structural changes due to fishing effects. J. Sea Res. 79, 50–59 (2013).Article 
ADS 

Google Scholar 
Mendonça, V. et al. What’s in a tide pool? Just as much food web network complexity as in large open ecosystems. PLoS ONE 13, e0200066 (2018).Article 

Google Scholar 
Mahon, R. & Mahon, S. D. Structure and resilience of a tidepool fish assemblage at Barbados. Environ. Biol. Fishes 41, 171–190 (1994).Article 

Google Scholar 
Willis, T. J. & Roberts, C. D. Recolonisation and recruitment of fishes to intertidal rockpools at Wellington, New Zealand. Environ. Biol. Fish. 47, 329–343 (1996).Article 

Google Scholar 
Ho, L.-T., Wang, S.-C., Shao, K.-T., Chen, I.-S. & Chen, H. A long-term monitoring dataset of fish assemblages in rocky tidepools on the northern coast of Taiwan. Sci. Data 7, 84 (2020).Article 

Google Scholar 
Chen, H. et al. Long-term monitoring dataset of fish assemblages impinged at nuclear power plants in northern Taiwan. Sci. Data 2, 150071 (2015).Article 
CAS 

Google Scholar 
Liao, Y. C., Chen, L. S., Shao, K. T. & Tu, Y. Y. Temporal changes in fish assemblage from the impingement data at the second nuclear power plant, northern Taiwan. J. Mar. Sci. Technol. 12, 411–417 (2004).Article 

Google Scholar 
Chen, H., Chen, C.-Y. & Shao, K.-T. Recovery and variation of the coastal fish community following a cold intrusion event in the Penghu Islands, Taiwan. Plos One 15, e0238550 (2020).Article 
CAS 

Google Scholar 
Heard, J. et al. Coastal development threatens Datan area supporting greatest fish diversity at Taoyuan Algal Reef, northwestern Taiwan. Aquat. Conserv.: Mar. Freshw. Ecosyst. 31, 590–604 (2021).Article 

Google Scholar 
Chen, H., Shao, K. T. & Kishino, H. Phylogenetic skew: an index of community diversity. Mol. Ecol. 24, 759–770 (2015).Article 

Google Scholar 
Chen, H., Shao, K.-T. & Kishino, H. Bayesian hierarchical ANOVA model of stochastic seasonality for Diodon holocanthus in northern Taiwan. J. Mar. Sci. Technol. 24, 303–310 (2016).
Google Scholar 
Chen, H., Chen, C.-Y. & Shao, K.-T. Time series dataset of fish assemblages near thermal discharges at nuclear power plants in northern Taiwan. Sci. Data 5, 180055 (2018).Article 

Google Scholar 
Chen, H., Chen, C.-Y., Shao, K.-T. & Gong, G.-C. Spatial and temporal variations in species diversity of fish assemblages near a sewage treatment plant in northern Taiwan. Fish. Sci. 85, 581–590 (2019).Article 
CAS 

Google Scholar 
Shao, K.-T., Chen, J.-P. & Wang, S.-C. in Proceedings of the 5th Indo-Pacific Fish Conference (eds. Seret, B. & Sire, J.-Y.) Biogeography and database of marine fishes in Taiwan waters (French Society of Ichthyology, 1999).Shao, K.-T. et al. A checklist of the fishes of southern Taiwan, northern South China Sea. Raffles Bull. Zool. Suppl. 19, 233–271 (2008).
Google Scholar 
Shen, S.-C. Chief Editor. Fishes of Taiwan. (Department of Zoology, National Taiwan University, 1993).Ho, L-T., Wang, S-C., Shao, K-T., Chen, I-S. & Chen, H. Long-term monitoring dataset of fish assemblages in rocky tidepools on the southern coast of Taiwan, figshare, https://doi.org/10.6084/m9.figshare.20747467.v1 (2022). More

We applied a modified version of a planning framework for CWR conservation25,26 which has been used by numerous countries of Europee.g.29,63,64, Americae.g.65, Africa30 and Asia66,67. We addressed the following main steps of the toolkit (see Spanish version49): (i) CWR checklist, i.e., creating a list of CWR taxa distributed in an area (Supplementary Data 1), (ii) CWR inventory, i.e., taxa selection and collation of ancillary data, including taxonomic data (Supplementary Data 2), (iii) taxa extinction risk assessment (Table 1, Supplementary Data 3), and (iv) a systematic conservation planning assessment, i.e., spatial analyses to assess conservation areas (Fig. 1). We only provide a brief description of steps i-iii, as these are thoroughly described in Goettsch et al.2. Here, we focus on the systematic conservation planning assessment, introducing an approach in order to identify conservation areas for CWR that account for genetic differentiation in a spatially explicit way, through the use of proxies of genetic differentiation (Fig. 1).During the process -framed under the project “Safeguarding Mesoamerican crop wild relatives” (https://www.darwininitiative.org.uk/project/23007/)- more than 100 experts from academic, governmental, and non-governmental organizations from El Salvador, Guatemala, Honduras, Mexico, the UK, and IUCN participated in six workshops, shared data, and provided fundamental knowledge and feedback at each project stage to ensure accurate, reliable and robust information for next steps. The checklist, inventory and risk assessment were collaboratively developed between partners of El Salvador, Guatemala, and Mexico (hereafter, Mesoamerica; Goettsch et al.2). The spatial analysis to identify areas for in situ and ex situ conservation of CWR was done independently by each country.To assess conservation areas of CWR in Mexico, we developed proxies of genetic differentiation that account for evolutionary processes by including historical and environmental drivers of genetic diversity (see the Methods section ‘Proxies of genetic differentiation’). In addition, we used criteria such as information on taxon-specific tolerance to human-modified habitats and IUCN extinction risk category. We applied a systematic conservation planning approach and performed spatial analysis using the software Zonation50. We compared different scenarios to represent genetic diversity of CWR based on potential species distribution models (SDM) and proxies of genetic differentiation.Study areaMesoamerica is a cultural region encompassing the territories of Belize, Guatemala, El Salvador, the southern part of Mexico and parts of Honduras, Nicaragua and Costa Ricasee 2. In this study, we also included the dry areas of northern Mexico that are part of Aridamerica68 and the Nearctic biogeographic realm69 to account for the full extent of the geographic range of many taxa included in the extinction risk assessment2.For the assessment of conservation areas, we focused on Mexico, which is one of the most biodiverse countries in the world70. The Mexican territory covers 80% of the landscapes of the region called Mesoamerica. Its high biological diversity is attributed to its geographic, topographic, climatic, geological and cultural characteristics, which, among other factors, shaped the distribution of an extraordinary variety of ecosystems and species with high levels of endemism and species turnover among different regions32,71,72,73. In particular, the high genetic variation within populations of landraces and CWR is the result of past and ongoing sociocultural processes occurring in a wide range of distinct environmental conditions74,75.(i) CWR checklist and (ii) CWR inventoryThe compiled CWR checklist included ~3000 species and subspecies of 92 genera and 45 families of plants that belong to the same genus of a crop cultivated in Mesoamerica, or wild plant collected for food or other uses in the region (Supplementary Data 1).The first set of criteria were established in preparation for the first stakeholder workshop. The following criteria were applied at the genus level to compile the CWR inventory: (1) occurrence of wild relatives of cultivated plants or crops that were domesticated in Mesoamerica; (2) existence of research groups working on taxa that could support the extinction risk assessment; and (3) relation to a crop of economic and nutritional importance at local, national and regional levels, or cultivars known to require genetic improvement.To narrow the list for the inventory and extinction risk assessment, similar criteria were agreed upon in the same workshop and applied at the species level: (1) native distribution in Mesoamerica, incl. Aridamerica; (2) related to a crop of economic or social importance based on production and nutritional value; (3) related to a taxon for which Mesoamerica is the center of origin or domestication; (4) constitutes part of the primary or secondary gene pool, and in some cases the tertiary gene pool76. The primary gene pool consists of wild plants of the same species as the crop and thus their mating produces strong fertile progeny. The secondary gene pool is composed of wild relatives distinct from cultivated species but closely related as to produce some fertile offspring (same taxonomic series or section in the absence of crossing and genetic diversity information, see the ‘taxon group’ concept proposed by Maxted and collaborators77, Supplementary Note 5). The tertiary gene pool (same subgenus in the taxon group concept) corresponds to CWR that are more distant relatives to the taxa of the primary gene pool, but can have important adaptive traits which can be used with specific breeding techniques. This provided a preliminary list of 514 CWR taxa related to avocado, cotton, amaranth, cocoa, squash, sweet potato, chayote, chili pepper, cempasuchil, bean, sunflower, maize, papaya, potato, vanilla, and yuca (Supplementary Data 2).The list had to be further reduced due to time and funding restrictions to include those genera which when added together would include no more than 250 taxa, and that the taxonomic groups could be comprehensively assessed and their taxa evaluated throughout their entire range. Thus, not all species in the group necessarily met the criteria previously mentioned. See the final Mesoamerican CWR inventory in Supplementary Data 3; see summary in Table 1.(iii) Taxa extinction risk assessmentFull methodological details and results of this section are described in Goettsch et al.2. Summarizing, during the process 224 taxa were evaluated according to the International Union for Conservation of Nature, IUCN, Red List Categories and Criteria78. The IUCN Red List is a critical indicator to identify species most vulnerable to extinction considering a set of criteria, i.e., species’ population trends, size, structure, and geographic ranges. A Red List workshop with the participation of 25 experts from different project partner institutions and IUCN specialists was organized to assess the extinction risk of taxa. The threat analysis included not only species, but subspecies and subpopulations (i.e. races) for some groups (Supplementary Data 3, see summary in Table 1).(iv) Systematic conservation planning assessmentTo undertake the following spatial analyses we focused on the dataset of 224 CWR described above, which is representative of the CWR of the main crops of Mesoamerica (10 genera, Table 1).Species distribution modelingTo compile occurrence records, hundreds of data sources were consulted, including published and personal databases of the project participantse.g.79,80,81,82, the Agrobiodiversity Atlas of Guatemala (https://www.ars.usda.gov/northeast-area/beltsville-md-barc/beltsville-agricultural-research-center/national-germplasm-resources-laboratory/docs/atlas-of-guatemalan-crop-wild-relatives), the Global Biodiversity Information Facility (GBIF, https://www.gbif.org/), and Mexico’s Biodiversity Information System (SNIB, http://snib.mx/).To generate potential species distribution models (SDM), we used more than 13,000 occurrence records (Supplementary Data 4), that were standardized and curated by experts to generate the range maps of taxa as part of the extinction risk assessment, which were published in IUCN Red List (https://www.iucn.org/news/species/202109/threats-crop-wild-relatives-compromising-food-security-and-livelihoods). Spatial resolution of the SDM was 1 km2. SDM were obtained for taxa with more than 20 unique occurrence data in a 1 km2 grid covering the study extent to reduce uncertainty when using smaller sample sizes83. We used 19 bioclimatic variables and other climatic variables, such as annual potential evapotranspiration, aridity index, annual radiation, slope, and altitude84,85,86. Climate data represents annual and seasonal patterns of climate between 1950 and 2000. Also, we used a variable that described the percentage of bare soil and cultivated areas87. Collinearity between variables was assessed with the ‘corselect’ function of the package fuzzySim version 1.088, using a value of 0.8 and the variance inflation factors as criteria to exclude highly correlated variables.We used MaxEnt version 3.3.1, a machine-learning algorithm that uses the maximum entropy principle to identify a target probability distribution, subject to a set of constraints related to the occurrence records and environmental data89,90. Model calibration area for each taxon included those ecoregions where the taxon has been recorded; we used the terrestrial ecoregions dataset69. We did this based on the calibration area or ‘M element’ of the BAM diagram that refers to areas that have been accessible to the taxon via dispersal over relevant periods of time91,92. We randomly sampled 10,000 background localities from the selected areas.To reduce model complexity without compromising model performance, we built several models by varying the feature classes (FC) and regularization multipliers (RM) (see refs. 93,94,95) using R 3.6.096 and ‘ENMeval’ version 0.3.0 package97. FC determines the flexibility of the modeled response to the predictor variables, while the RM penalizes model complexity93. Occurrence records were randomly divided into 70% for model selection, and 30% of data was withheld for model validation. ENMeval carries out an internal partition of localities to test each combination of settings. Therefore, we selected the random k-fold method to divide localities into four bins. We build models with six FC combinations and varied RM values ranging from 0.5 to 4.0 in 0.5 increments. Optimal models were selected using Akaike’s Information Criterion corrected for small sample sizes (⍙AICc = 0). This method penalizes overly complex models and helps to choose those with an optimal number of parameters. However, it has been shown that the number of model parameters may not correctly estimate degrees of freedom98, and that model selection should not be selected solely with one measure99. Thus, we used 30% of the withheld data to test the area under the curve (AUC) of the receiver operating characteristic, and the omission error under a 10 percentile training threshold.We used the ten percentile or minimum training presence threshold to obtain binary maps of the presence and absence of suitable areas for species distribution. We asked experts of each taxonomic group who were also involved in the extinction risk assessment to select one of these two options and to indicate possible overestimated areas, which were then eliminated case by case using the information of Mexican ecoregions100 and watersheds101. Eight models were binarized with the minimum training presence threshold; for the other models we used the 10 percentile threshold. See MaxEnt performance and significance of SDM at Supplementary Data 5. AUC values ranged from 0 to 1; 0.5 indicated a model performance not better than random, while values closer to 1 indicated a better model performance; here we used SDM showing AUC values higher than 0.7. For Phaseolus and Zea, we used SDM that were previously generated by Delgado-Salinas et al.102, and Sánchez González et al.103, respectively. SDM for 116 taxa were validated by experts of each taxonomic group. See references and download links at Supplementary Data 6.For the conservation planning analysis of Mexico, we clipped the models to the Mexican territory, and trimmed the continuous SDM using the binary SDM to keep pixel values of areas with elevated probability of taxa presence. For taxa without SDM, we included the occurrence records of these taxa in the spatial analysis by using the information on observation location, i.e., coordinates (see Supplementary Data 3). This is done by enabling the function ‘species of special interest’ (SSI). See further details in the method section ‘Final conservation analysis’.Proxies of genetic differentiationTo identify proxies of genetic differentiation in an explicit, efficient, and repeatable way, we included environmental and historical drivers of genetic diversity. For this, we first divided Mexico into 27 Holdridge life zones (Supplementary Fig. 2, Supplementary Data 8), which we then subdivided according to phylogeographic studies that have found genetic differentiation among populations of several taxa (see division of each life zone into proxies in Supplementary Fig. 4; Supplementary Fig. 3 provides a general geographical overview of Mexico and main geographic references mentioned in Supplementary Fig. 4). The literature review was done searching for the words “phylogeography” and one of the following: (i) name of the Mexican biogeographic zones, (ii) “Mexico” + an ecosystem name (e.g. “Mexico” “rainforest”) or (iii) “Mexico” + lowlands/highlands. See list of references used in this study in Supplementary Data 9.In addition, we manually reviewed the citations to the most cited papers of the previous search. Reviews and meta-analyses were also included, although we excluded studies performed in CWR to show that our approach can be used without prior information on this group. As more studies on such taxa become available, they can be used to fine-tune the proxies of genetic differentiation. We focused on terrestrial species including plants, animals, and fungi (Supplementary Data 10) except to subdivide a life zone covering the coasts of the California Peninsula, where we could not find studies on terrestrial taxa so we included studies on fish species (see Supplementary Fig. 4).Since most of the life zones cover large territories, and complete phylogeographic congruence among different taxa is uncommon, we targeted to represent general trends that would likely occur across diverse species, instead of trying to represent fine idiosyncratic patterns of genetic differentiation. For instance, although distribution ranges of highland taxa shifted during the Pleistocene climate fluctuations, in general populations persisted (glacial-interglacial periods) within the main mountain ranges, while lowland populations were ephemeral (only glacial periods). So, gene flow among mountain ranges was more limited than within them. As a result, genetic differentiation among mountain ranges of different biogeographic provinces has been widely documented32, so we used this general pattern to subdivide the life zones that occur in highlands. These types of patterns are particularly relevant for a country like Mexico, due to its complex topography, tropical latitude, and geographic features of different ages, which promote population differentiation among the Mexican main geographic features. To translate the phylogeographic information into a spatial context, we used biogeographic regions, basins, topographic or edaphic data to split the life zones into different subzones using the best fitting cartography to represent the phylogeographic patterns (Supplementary Fig. 4).We obtained 102 proxies of genetic differentiation for Mexico (Supplementary Fig. 5). We validated our findings by using available genomic data of an empirical study of a wild relative of maize, the teosinte Zea mays subsp. parviglumis, which was not included in the literature review in order to test the usefulness of our approach regarding the lack of genetic data. The dataset includes ca. 1800 occurrence records and ca. 30,000 SNPs48. Sampling localities were not used for distribution modeling. Admixture groups per population were estimated for K1 to 60. According to the population analysis, Z. mays subsp. parviglumis is structured in 13 genetic clusters along a longitudinal gradient (Fig. 3a–c). We used the K = 13 for plotting based on the Cross-Validation error. The proportion of each genetic cluster was estimated by sampling locality and plotted using pie charts over the map (Supplementary Fig. 6). Then, using the data layer of the SDM subdivided by proxies of genetic differentiation, we extracted which was the proxy most frequent in a 5 km buffer for each sampling locality. The Admixture plot was ordered by all genetic clusters and subdivided by the proxy of genetic differentiation most frequent for each locality. In addition, we calculated a principal component analysis (PCA) and projected into a score plot the first three components. Individual samples were colored by the proxies where they fell in the 5 km buffer (Fig. 3c). To compare how genetic variation was represented by the different scenarios we plotted the proportion of the area of each proxy as given by the potential SDM according to two different scenarios (only considering SDM; combining SDM*PGD) considering 20% of Mexico’s terrestrial area (Fig. 3d). Analyses were run in R version 3.5.196 using the R packages pcadapt version 4.3.3104, ggplot2 version 2_3.3.3105, readr version 1.4.0106, gridExtra version 2.3107, ggnewscale version 0.4.5108, scatterpie version 0.1.5109, pophelper version 2.3.1110, raster version 3.4-5111, rgdal version 1.4-8112, rgl version 0.107.10113, and sp version 1.4-4114,115.Habitat preferenceWe considered habitat preference to refine the presence of CWR in the planning process; thus minimizing commission errors and highlighting areas that more probably contain taxa116. For each taxon, experts assessed its habitat preference (1: high preference; 0.5: low preference; 0.1: no preference) according to the following categories: (i) well-conserved vegetation (i.e. primary vegetation), (ii) human-impacted vegetation (i.e. secondary vegetation), (iii) less intensive rainfed and moisture agriculture, (iv) intensive rainfed and moisture agriculture, (v) irrigated agriculture, (vi) induced and cultivated grasslands and forests, and vii) urban areas (Supplementary Data 11). To spatially delimit these classes, we used the land use cover and vegetation map for Mexico117, and assessed seven main categories of land cover by grouping the map legend (Supplementary Fig. 9). To differentiate between less intensive and intensive cultivated areas, we followed Bellon et al.56, who associated the presence of native maize varieties of Mexico to occur in municipalities with average yields of less than or equal to 3 t ha-1 using agricultural production data from 2010 from the Information System of Agrifood and Fisheries (SIAP), and selected the municipalities with the established average maize yield. We combined the municipality layer with the land cover map to differentiate areas of high and low agricultural intensity. To generate taxon-specific habitat layers, we associated the habitat preference classes established by experts to the land cover map aggregated into seven major land cover categories, using R 3.6.096 and the following packages: raster version 3.4-5111 and rgdal version 1.4-8112. We obtained habitat maps for 116 taxa with SDM.Preliminary analysisWe generated five preliminary scenarios to explore different approaches to include conservation features for maximizing the representation of intraspecific diversity as given by taxa and proxies of genetic differentiation, i.e., representation of proxies within a taxa range (Supplementary Fig. 7): (i) “SDM” scenario, included 116 SDM, which we used as base scenario to examine the representation of taxa and proxies of genetic variability (n = 116); (ii) “SDM + LZ” scenario, included 116 SDM and 27 layers representing Holdridge life zones to consider environmental variation (n = 143); (iii) “SDM + PGD” scenario, included 116 SDM and 102 layers representing each proxy of genetic differentiation individually (n = 218); (iv) “SDM*PGD” scenario, included 5004 input layers representing the intersection of SDM and PGD (n = 5004; combining 116 SDM with 102 proxies resulted in 11,832 layers, but as some of the intersections produced empty outputs given the extension of SDM that do not cover all Mexico, for further analysis we used 5004 input layers with value data. To subdivide the layers, we used ArcGIS version 10.2.2118; to filter the layers, we used R 3.5.196.); (v) “SDM and PGD as ADMU” scenario, included 116 SDM as the main conservation features, while integrating one single layer of proxies of genetic differentiation to consider each of them as planning units by using the ‘Administrative units’ function. Analysis was done in Zonation50,119.We compared the results by assessing 20% of Mexico’s terrestrial area (Fig. 5b) to perform statistical analysis in R 3.5.196 using the following packages: purrr version 0.3.4120, ‘dplyr’ version 1.0.2121, ‘ggplot2’ version 2_3.3.3105, ‘raster’ version 3.4-5111, ‘scales’ version 1.2.0122, ‘sp’ version 1.4-4114,115, ‘tidyr’ version 1.0.2123, and ‘vegan’ version 2.6-2124. The area threshold was established based on Aichi target 11 and on comparisons of performance curves to efficiently represent taxa ranges delimited by SDM and proxies of genetic differentiation (Fig. 6). As using SDM combined with proxies of genetic differentiation showed the highest representation of genetic diversity (“SDM*PGD” scenario), we used this approach for the final analyses.Final conservation analysisWe identified areas of high conservation value for CWR in Mexico by using the software Zonation version 4.050,119, a systematic conservation planning tool that allows optimizing representation of species, taxa, or other conservation features, e.g., proxies of genetic differentiation, in a given study area. The program hierarchically ranks areas by removing cells of low conservation value, as given, for example, by a reduced number of taxa or occurrence of low weighted features, while considering multiple criteria such as the weighting of taxa and habitat preference of taxa. We applied the core-area zonation removal rule (CAZ) to maximize the representation of all conservation features in a minimal possible area51. Zonation generates two main outputs: (a) a hierarchical landscape priority rank map, that allows decision makers establishing different area thresholds to highlight areas of conservation interest; and (b) a representation curve showing species or conservation features range distribution in a given area. The curve also allows identifying how much area is needed to cover a certain taxon range or the distribution of a feature of conservation interest.For the conservation scenarios, we integrated the following inputs in the Zonation software: (1) 5,004 layers, i.e., SDM intersected with proxies of genetic differentiation (as described by “SDM*PGD” scenario, Fig. 4), (2) occurrence records of 98 taxa; only for those taxa without SDM, see Supplementary Data 3), (3) taxa specific habitat layers (according to Supplementary Data 11 and Supplementary Fig. 9), and (4) IUCN threat category (Supplementary Data 3) as an additional parameter to weight taxa differently to consider their vulnerability to extinction, see details below. See Zonation configuration at Supplementary Note 6.Data from different sources can be mixed in the same analysis, which is useful to not lose or omit information of any taxa of interest in the assessment. Here, we included information of a total of 214 taxa (see Supplementary Data 3). Distribution data of 116 taxa were represented by 5004 layers that resulted from combining 116 SDM and 102 PGD. This approach showed the highest proportion of area of taxa ranges (on average 41%) and highest representation of PGD within the area of each taxon (on average 76%; Fig. 4; see description in the main text). For some taxa, e.g. Cucurbita pepo, Physalis cinerascens, and Zea mays information on its distribution was assessed at subspecies level rather than at species level, explaining the difference in numbers of CWR taxa.In addition, we included occurrence data of 98 taxa without SDM to prevent missing important areas of taxa known distribution that are important to conserve (see Supplementary Data 3). We enabled the function ‘species of special interest’ (SSI) of Zonation, and included a SSI feature list file, listing the taxon names, as well as taxon-specific coordinate file for each of the 98 taxa that have been reviewed by the experts of each group. The spatial reference system was World Mercator projection. Occurrence data and SDM are treated similarly in the Zonation analysis, i.e., cells where taxa occur will be retained in the solution as long as possible to maximize its representation in the solution.We assigned weights to the 116 taxa with SDM by using IUCN threat categories (according to Supplementary Data 3), giving highest values to taxa with highest risk of extinction that urgently need management actions to further avoid genetic erosion. By including conservation feature weights, Zonation estimates the conservation value of a cell not only based on the presences of a taxa and their distribution range, but also on the weight. A high weight indicates a high conservation value of cells where these taxa are distributed. As there is no rule for weight setting, we assigned values between 1 and 0 regardless of taxa distribution ranges, which is automatically considered in the Zonation algorithm to guarantee the representation of locations where limited-range distributed taxa occur within the most valuable conservation area. Thus, weights were assigned as follows: Critically endangered, CR: 1; Endangered, EN: 1; Vulnerable, VU: 0.8; Near threatened, NT: 0.5; Data deficient, DD: 0.3; Least concern, LC: 0.2 Not evaluated, NE: 0.1. SSI taxa were all weighted similarly with 1 in order to represent the 98 SSI taxa and their occurrences in the top fraction of the most valuable conservation area, as these areas could be considered as ‘irreplaceable’ in terms of conservation. The conservation of these taxa that are only known in a few locations is crucial to maintain their populations. Information on weights for taxa with and without SDM is included in the file that lists the 5004 conservation features and the SSI file, respectively.To include the information on habitat, we included 116 habitat maps which guide the selection of cells to areas where its presence is more probable (see the Methods section: “Habitat preference”). This option can only be used for taxa represented by a raster layer, and is not available for SSI taxa included via occurrence records. By enabling the “landscape condition” option of Zonation, each habitat map is linked to a specific conservation feature layer. Areas with unfavorable habitats will quickly be masked out during the selection of cells in order to obtain a solution that favors conservation areas within areas of preferred habitat.We generated three final scenarios to identify conservation areas for (a) all taxa, (b) taxa exclusively distributing in natural vegetation, and (c) taxa associated with a wider range of habitats such as natural vegetation, agricultural and urban areas. The Zonation configuration remained similar among the three scenarios. When taxa were not included in a given scenario, we assigned a value weight of 0. This excluded the feature to be considered for the hierarchical prioritization of the landscape, but still allowed to evaluate the taxa during post-processing.To evaluate the spatial results (Supplementary Fig. 11), we analyzed performance curves to represent proxies of genetic differentiation within each taxon range (Supplementary Fig. 12). Also, we considered the most valuable 20% area of Mexico to calculate the coincidence of the three scenarios (Supplementary Fig. 13), and the overlap with federal protected areas125 and indigenous regions126,127 (Supplementary Fig. 14), and land cover data used in the analyses (Supplementary Figs. 9, 15).We discussed the proposed methodological framework, input layer and criteria during a fourth workshop in Mexico. It is worth mentioning that we ran several analyses including additional layers, such as areas where indigenous communities live that promote the presence of CWR in the landscape6. However, as the output indicated no evident difference by including this information, final analyses did not consider these data. We neither included protected areas nor tried to expand on the current 12% protected area system, because most management plans do not specifically address CWR management (but see the management program of the Protected Area of ‘Sierra de Manantlán’128), and thus generally do not adequately plan for wild and native genetic resources129. We also discussed different approaches to consider connectivity for taxa, habitats and proxies of genetic differentiation in the Zonation processing. Still, we finally decided to run the analysis without particularly accounting for connectivity as we had no taxa-specific information on dispersal abilities or possible effects of fragmentation, and we did not want to lose efficiency of the solution to represent taxa by or include lower-quality habitats by forcing the solution to an aggregation of pixels.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

Blatchford, M. L., Mannaerts, C. M., Zeng, Y., Nouri, H. & Karimi, P. Status of accuracy in remotely sensed and in-situ agricultural water productivity estimates: A review. Remote Sensing of Environment 234, 111413, https://doi.org/10.1016/j.rse.2019.111413 (2019).Article 
ADS 

Google Scholar 
Geerts, S. & Raes, D. Deficit irrigation as an on-farm strategy to maximize crop water productivity in dry areas. Agricultural Water Management 96, 1275–1284, https://doi.org/10.1016/j.agwat.2009.04.009 (2009).Article 

Google Scholar 
Hellegers, P., Soppe, R., Perry, C. & Bastiaanssen, W. Combining remote sensing and economic analysis to support decisions that affect water productivity. Irrigation Science 27, 243–251, https://doi.org/10.1007/s00271-008-0139-7 (2009).Article 

Google Scholar 
Bastiaanssen, W. G. M. & Steduto, P. The water productivity score (WPS) at global and regional level: Methodology and first results from remote sensing measurements of wheat, rice and maize. The Science of the total environment 575, https://doi.org/10.1016/j.scitotenv.2016.09.032 (2017).Seneviratne, S. I. et al. Investigating soil moisture–climate interactions in a changing climate: A review. Earth Science Reviews 99, https://doi.org/10.1016/j.earscirev.2010.02.004 (2010).Hu, X., Shi, L., Lin, L. & Zha, Y. Nonlinear boundaries of land surface temperature–vegetation index space to estimate water deficit index and evaporation fraction. Agricultural and Forest Meteorology 279, https://doi.org/10.1016/j.agrformet.2019.107736 (2019).Bowen, I. S. The Ratio of Heat Losses by Conduction and by Evaporation from any Water Surface. Physical Review 27, 779–787, https://doi.org/10.1103/PhysRev.27.779 (1926).Article 
ADS 
CAS 
MATH 

Google Scholar 
Penman, H. L. Natural evaporation from open water, hare soil and grass. Proceedings of the Royal Society of London. Series A, Mathematical and physical sciences 193, https://doi.org/10.1098/rspa.1948.0037 (1948).Monteith, J. L. Evaporation and environment. The stage and movement of water in living organisms. Symp.soc.exp.biol.the Company of Biologists (1965).Wang, K. & Dickinson, R. E. A review of global terrestrial evapotranspiration: Observation, modeling, climatology, and climatic variability. Reviews of Geophysics 50, https://doi.org/10.1029/2011RG000373 (2012).Bastiaanssen, W. G. et al. A remote sensing surface energy balance algorithm for land (SEBAL) Part 1: Fomulation. Journal of hydrology 212, 213–229, https://doi.org/10.1016/S0022-1694(98)00253-4 (1998).Article 
ADS 

Google Scholar 
Bastiaanssen, W. G. M. et al. A remote sensing surface energy balance algorithm for land (SEBAL) Part 2. Validation. Journal of Hydrology 212, https://doi.org/10.1016/S0022-1694(98)00254-6 (1998).Su, Z. The Surface Energy Balance System (SEBS) for estimation of turbulent heat fluxes. Hydrology and Earth System Science 6, 85–99, https://doi.org/10.5194/hess-6-85-2002 (2002).Article 
ADS 

Google Scholar 
Norman, J. M., Kustas, W. P. & Humes, K. S. Source approach for estimating soil and vegetation energy fluxes in observations of directional radiometric surface temperature. Agricultural and Forest Meteorology 77, https://doi.org/10.1016/0168-1923(95)02265-y (1995).Mu, Q., Heinsch, F. A., Zhao, M. & Running, S. W. Development of a global evapotranspiration algorithm based on MODIS and global meteorology data. Remote Sensing of Environment 111, https://doi.org/10.1016/j.rse.2007.04.015 (2007).Mu, Q., Zhao, M. & Running, S. W. Improvements to a MODIS global terrestrial evapotranspiration algorithm. Remote Sensing of Environment 115, 1781–1800, https://doi.org/10.1016/j.rse.2011.02.019 (2011).Article 
ADS 

Google Scholar 
Fisher, J. B., Tu, K. P. & Baldocchi, D. D. Global estimates of the land–atmosphere water flux based on monthly AVHRR and ISLSCP-II data, validated at 16 FLUXNET sites. Remote Sensing of Environment 112, 901–919, https://doi.org/10.1016/j.rse.2007.06.025 (2008).Article 
ADS 

Google Scholar 
Kim, H. W., Hwang, K., Mu, Q., Lee, S. O. & Choi, M. Validation of MODIS 16 global terrestrial evapotranspiration products in various climates and land cover types in Asia. KSCE Journal of Civil Engineering 16, https://doi.org/10.1007/s12205-012-0006-1 (2012).Velpuri, N. M., Senay, G. B., Singh, R. K., Bohms, S. & Verdin, J. P. A comprehensive evaluation of two MODIS evapotranspiration products over the conterminous United States: Using point and gridded FLUXNET and water balance ET. Remote Sensing of Environment 139, https://doi.org/10.1016/j.rse.2013.07.013 (2013).Jin, X. et al. Estimation of water productivity in winter wheat using the AquaCrop model with field hyperspectral data. Precision Agriculture 19, 1–17, https://doi.org/10.1007/s11119-016-9469-2 (2016).Article 

Google Scholar 
Felix, R., Clement, A., Igor, S. & Oscar, R. Using Low Resolution Satellite Imagery for Yield Prediction and Yield Anomaly Detection. Remote Sensing 5, 1704–1733, https://doi.org/10.3390/rs5041704 (2013).Article 

Google Scholar 
Lu, Y. et al. Assimilation of soil moisture and canopy cover data improves maize simulation using an under-calibrated crop model. Agricultural Water Management 252, https://doi.org/10.1016/j.agwat.2021.106884 (2021).Jin, X., Kumar, L., Li, Z., Feng, H. & Wang, J. A review of data assimilation of remote sensing and crop models. European Journal of Agronomy 92, https://doi.org/10.1016/j.eja.2017.11.002 (2018).Weiss, M., Jacob, F. & Duveiller, G. Remote sensing for agricultural applications: A meta-review. Remote Sensing of Environment 236, https://doi.org/10.1016/j.rse.2019.111402 (2019).Jin, X. et al. Winter wheat yield estimation based on multi-source medium resolution optical and radar imaging data and the AquaCrop model using the particle swarm optimization algorithm. ISPRS Journal of Photogrammetry and Remote Sensing 126, 24–37 (2017).Article 
ADS 

Google Scholar 
Tao, F., Rötter, R. P., Palosuo, T., Díaz-Ambrona, C. G. H. & Schulman, A. H. Contribution of crop model structure, parameters and climate projections to uncertainty in climate change impact assessments. Global Change Biology 24, https://doi.org/10.1111/gcb.14019 (2017).Jin, X. et al. A review of data assimilation of remote sensing and crop models. European Journal of Agronomy 92, 141–152, https://doi.org/10.1016/j.eja.2017.11.002 (2018).Article 

Google Scholar 
Anikó, K. et al. Statistical modelling of crop yield in Central Europe using climate data and remote sensing vegetation indices. Agricultural and Forest Meteorology 260-261, 300–320, https://doi.org/10.1016/j.agrformet.2018.06.009 (2018).Article 

Google Scholar 
Wang, Y., Zhang, Z., Feng, L., Du, Q. & Runge, T. Combining Multi-Source Data and Machine Learning Approaches to Predict Winter Wheat Yield in the Conterminous United States. Remote Sensing 12, 1232, https://doi.org/10.3390/rs12081232 (2020).Article 
ADS 

Google Scholar 
Franz, T. E. et al. The role of topography, soil, and remotely sensed vegetation condition towards predicting crop yield. Field Crops Research 252, https://doi.org/10.1016/j.fcr.2020.107788 (2020).Noland, R. L. et al. Estimating alfalfa yield and nutritive value using remote sensing and air temperature. Field Crops Research 222, 189–196, https://doi.org/10.1016/j.fcr.2018.01.017 (2018).Article 

Google Scholar 
Cao, J., Zhang, Z., Luo, Y., Zhang, L. & Tao, F. Wheat yield predictions at a county and field scale with deep learning, machine learning, and google earth engine. European Journal of Agronomy, 126204, https://doi.org/10.1016/j.eja.2020.126204 (2021).Jacinta, H. & Kerrie, M. Statistical Machine Learning Methods and Remote Sensing for Sustainable Development Goals: A Review. Remote Sensing 10, 1365, https://doi.org/10.3390/rs10091365 (2018).Article 

Google Scholar 
Jin, X., Liu, S., Baret, F., Hemerlé, M. & Comar, A. Estimates of plant density of wheat crops at emergence from very low altitude UAV imagery. Remote Sensing of Environment 198, 105–114, https://doi.org/10.1016/j.rse.2017.06.007 (2017).Article 
ADS 

Google Scholar 
Maimaitijiang, M. et al. Soybean yield prediction from UAV using multimodal data fusion and deep learning. Remote Sensing of Environment 237, 111599, https://doi.org/10.1016/j.rse.2019.111599 (2020).Article 
ADS 

Google Scholar 
Hossein, A., Mohsen, A., Davoud, A., Salehi, S. H. & Soheil, R. Machine Learning Regression Techniques for the Silage Maize Yield Prediction Using Time-Series Images of Landsat 8 OLI. IEEE Journal of Selected Topics in Applied Earth Observations Remote Sensing PP, 1–15, https://doi.org/10.1109/JSTARS.2018.2823361 (2018).Johansen, K. et al. Predicting Biomass and Yield in a Tomato Phenotyping Experiment Using UAV Imagery and Random Forest. Frontiers in Artificial Intelligence 3, 28, https://doi.org/10.3389/frai.2020.00028 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
Zhang, L., Ding, X., Shen, Y., Wang, Z. & Wang, X. Spatial Heterogeneity and Influencing Factors of Agricultural Water Use Efficiency in China. Resources and Environment in the Yangtze Basin 28, https://doi.org/10.11870/cjlyzyyhj201904008 (2019).Cheng, M. et al. Satellite time series data reveal interannual and seasonal spatiotemporal evapotranspiration patterns in China in response to effect factors. Agric. Water Manage. 255, https://doi.org/10.1016/j.agwat.2021.107046 (2021).Zhou, L. Comprehensive agricultural regionalization in China. (Agricultural Press of China, 1985).Luo, Y., Zhang, Z., Chen, Y., Li, Z. & Tao, F. ChinaCropPhen1km: A high-resolution crop phenological dataset for three staple crops in China during 2000-2015 based on LAI products. Figshare https://doi.org/10.6084/m9.figshare.8313530.v6 (2019).Luo, Y., Zhang, Z., Chen, Y., Li, Z. & Tao, F. ChinaCropPhen1km: a high-resolution crop phenological dataset for three staple crops in China during 2000–2015 based on leaf area index (LAI) products. Earth System Science Data 12, 197–214, https://doi.org/10.5194/essd-12-197-2020 (2020).Article 
ADS 

Google Scholar 
Song, D. Second China Soil Survey. (Chinese Science Press, 1979).Zhang, T., Yang, X., Wang, H., Li, Y. & Ye, Q. Climatic and technological ceilings for Chinese rice stagnation based on yield gaps and yield trend pattern analysis. Global Change Biology 20, 1289–1298, https://doi.org/10.1111/gcb.12428 (2014).Article 
ADS 
PubMed 

Google Scholar 
Chen, Y., Zhang, Z. & Tao, F. Improving regional winter wheat yield estimation through assimilation of phenology and leaf area index from remote sensing data. European Journal of Agronomy 101, 163–173, https://doi.org/10.1016/j.eja.2018.09.006 (2018).Article 

Google Scholar 
Cheng, M. et al. Combining multi-indicators with machine-learning algorithms for maize yield early prediction at the county-level in China. Agricultural and Forest Meteorology 323, https://doi.org/10.1016/j.agrformet.2022.109057 (2022).Amir, J. & Sinclair, T. A model of the temperature and solar-radiation effects on spring wheat growth and yield. Field Crops Research 28, 47–58, https://doi.org/10.1016/0378-4290(91)90073-5 (1991).Article 

Google Scholar 
Prince, S. D., Haskett, J., Steininger, M. & Wright, S. R. Net Primary Production of U.S. Midwest Croplands from Agricultural Harvest Yield Data. Ecological Applications 11, 1194–1205, https://doi.org/10.1890/1051-0761(2001)011[1194:NPPOUS]2.0.CO;2 (2001).Article 

Google Scholar 
Gilardelli, C. et al. Downscaling rice yield simulation at sub-field scale using remotely sensed LAI data. European journal of agronomy 103, 108–116, https://doi.org/10.1016/j.eja.2018.12.003 (2019).Article 

Google Scholar 
Shakoor, R., Hassan, M. Y., Raheem, A. & Wu, Y.-K. Wake effect modeling: A review of wind farm layout optimization using Jensen׳ s model. Renewable and Sustainable Energy Reviews 58, 1048–1059, https://doi.org/10.1016/j.rser.2015.12.229 (2016).Article 

Google Scholar 
Breiman, L. Random Forests. Machine Learning https://doi.org/10.1023/A:1010933404324 (2001).Article 
MATH 

Google Scholar 
Li, L. et al. Crop yield forecasting and associated optimum lead time analysis based on multi-source environmental data across China. Agricultural and Forest Meteorology 308–309, https://doi.org/10.1016/j.agrformet.2021.108558 (2021).Wang, L. A., Zhou, X., Zhu, X., Dong, Z. & Guo, W. Estimation of biomass in wheat using random forest regression algorithm and remote sensing data. The Crop Journal 4, 212–219, https://doi.org/10.1016/j.cj.2016.01.008 (2016).Article 

Google Scholar 
Feng, P. et al. Dynamic wheat yield forecasts are improved by a hybrid approach using a biophysical model and machine learning technique. Agricultural and Forest Meteorology 285-286, 107922, https://doi.org/10.1016/j.agrformet.2020.107922 (2020).Article 
ADS 

Google Scholar 
Lu, F., Sun, Y. & Hou, F. Using UAV Visible Images to Estimate the Soil Moisture of Steppe. Water 12, 2334, https://doi.org/10.3390/w12092334 (2020).Article 
CAS 

Google Scholar 
Wang, S. et al. High spatial resolution monitoring land surface energy, water and CO2 fluxes from an Unmanned Aerial System. Remote Sensing of Environment 229, 14–31, https://doi.org/10.1016/j.rse.2019.03.040 (2019).Article 
ADS 

Google Scholar 
Chen, Y. et al. Comparison of satellite-based evapotranspiration models over terrestrial ecosystems in China. Remote Sensing of Environment 140, 279–293, https://doi.org/10.1016/j.rse.2013.08.045 (2014).Article 
ADS 

Google Scholar 
Peralta, N., Assefa, Y., Du, J., Barden, C. & Ciampitti, I. Mid-Season High-Resolution Satellite Imagery for Forecasting Site-Specific Corn Yield. Remote Sensing 8, 848, https://doi.org/10.3390/rs8100848 (2016).Article 
ADS 

Google Scholar 
Russello, H. Convolutional neural networks for crop yield prediction using satellite images. IBM Center for Advanced Studies (2018).You, J., Li, X., Low, M., Lobell, D. & Ermon, S. in Proceedings of the AAAI Conference on Artificial Intelligence.Moran, P. A. Notes on continuous stochastic phenomena. Biometrika 37, 17–23 (1950).Article 
MathSciNet 
CAS 

Google Scholar 
Imran, M., Stein, A. & Zurita-Milla, R. Using geographically weighted regression kriging for crop yield mapping in West Africa. International Journal of Geographical Information Systems 29, 234–257, https://doi.org/10.1080/13658816.2014.959522 (2015).Article 

Google Scholar 
Harries, K. Extreme spatial variations in crime density in Baltimore County, MD. Geoforum 37, 404–416, https://doi.org/10.1016/j.geoforum.2005.09.004 (2006).Article 

Google Scholar 
Ghulam, A. et al. Remote Sensing Based Spatial Statistics to Document Tropical Rainforest Transition Pathways. Remote Sensing 7, 6257–6279, https://doi.org/10.3390/rs70506257 (2015).Article 
ADS 

Google Scholar 
Maimaitijiang, M., Ghulam, A., Sandoval, J. S. O. & Maimaitiyiming, M. Drivers of land cover and land use changes in St. Louis metropolitan area over the past 40 years characterized by remote sensing and census population data. International Journal of Applied Earth Observation Geoinformation 35, 161–174, https://doi.org/10.1016/j.jag.2014.08.020 (2015).Article 
ADS 

Google Scholar 
Cheng, M. Long time series (2001-2015) high-resolution crop yield and water productivity dataset of China, Zenodo, https://doi.org/10.5281/zenodo.5121842 (2021).Martens, B., Miralles, D. G., Lievens, H., Schalie, R. D. & Verhoest, N. GLEAM v3: Satellite-based land evaporation and root-zone soil moisture. Geoscientific Model Development 10, https://doi.org/10.5194/gmd-10-1903-2017 (2017).Wang, W., Cui, W., Wang, X. & Chen, X. Evaluation of GLDAS-1 and GLDAS-2 forcing data and Noah model simulations over China at the monthly scale. Journal of Hydrometeorology 17, 2815–2833, https://doi.org/10.1175/JHM-D-15-0191.1 (2016).Article 
ADS 

Google Scholar 
Chen, X. et al. Development of a 10-year (2001–2010) 0.1° data set of land-surface energy balance for mainland China. Atmospheric Chemistry and Physics 14, 14471–14518, https://doi.org/10.5194/acp-14-13097-2014 (2014).Article 
ADS 
CAS 

Google Scholar 
Ramoelo, A. et al. Validation of Global Evapotranspiration Product (MOD16) using Flux Tower Data in the African Savanna, South Africa. Remote Sensing 6, https://doi.org/10.3390/rs6087406 (2014).Yang, X., Yong, B., Ren, L., Zhang, Y. & Long, D. Multi-scale validation of GLEAM evapotranspiration products over China via ChinaFLUX ET measurements. International Journal of Remote Sensing https://doi.org/10.1080/01431161.2017.1346400 (2017).Article 

Google Scholar 
Hu, G., Jia, L. & Menenti, M. Comparison of MOD16 and LSA-SAF MSG evapotranspiration products over Europe for 2011. Remote Sensing of Environment 156, 510–526, https://doi.org/10.1016/j.rse.2014.10.017 (2015).Article 
ADS 

Google Scholar 
Khan, M. S., Liaqat, U. W., Baik, J. & Choi, M. Stand-alone uncertainty characterization of GLEAM, GLDAS and MOD16 evapotranspiration products using an extended triple collocation approach. Agricultural and Forest Meteorology 252, 256–268, https://doi.org/10.1016/j.agrformet.2018.01.022 (2018).Article 
ADS 

Google Scholar 
Glenn, E. P. et al. Scaling sap flux measurements of grazed and ungrazed shrub communities with fine and coarse-resolution remote sensing. Ecohydrology 1, 316–329, https://doi.org/10.1002/eco.19 (2008).Article 

Google Scholar 
Gamon, J. A. Reviews and Syntheses: optical sampling of the flux tower footprint. Biogeosciences 12, 4509–4523, https://doi.org/10.5194/bg-12-4509-2015 (2015).Article 
ADS 

Google Scholar 
Cai, Y. et al. Integrating satellite and climate data to predict wheat yield in Australia using machine learning approaches. Agricultural and Forest Meteorology 274, 144–159, https://doi.org/10.1016/j.agrformet.2019.03.010 (2019).Article 
ADS 

Google Scholar 
Chen, X. et al. Prediction of Maize Yield at the City Level in China Using Multi-Source Data. Remote Sensing 13, https://doi.org/10.3390/rs13010146 (2021).Guo, Y. et al. Integrated phenology and climate in rice yields prediction using machine learning methods. Ecological Indicators 120, 106935, https://doi.org/10.1016/j.ecolind.2020.106935 (2021).Article 

Google Scholar 
Yuan, W. et al. Estimating crop yield using a satellite-based light use efficiency model. Ecological Indicators 60, 702–709, https://doi.org/10.1016/j.ecolind.2015.08.013 (2016).Article 

Google Scholar 
Anandhi, A. Growing degree days – Ecosystem indicator for changing diurnal temperatures and their impact on corn growth stages in Kansas. Ecological Indicators 61, 149–158, https://doi.org/10.1016/j.ecolind.2015.08.023 (2016).Article 

Google Scholar 
Wart, J. V. Estimating Crop Yield Potential At National Scales. Field Crops Research 143, 34–43, https://doi.org/10.1016/j.fcr.2012.11.018 (2013).Article 

Google Scholar 
Kang, Y. S. et al. Yield prediction and validation of onion (Allium cepa L.) using key variables in narrowband hyperspectral imagery and effective accumulated temperature. Computers and Electronics in Agriculture 178, https://doi.org/10.1016/j.compag.2020.105667 (2020).Long, D., Singh, V. P. & Li, Z.-L. How sensitive is SEBAL to changes in input variables, domain size and satellite sensor? Journal of Geophysical Research: Atmospheres 116, https://doi.org/10.1029/2011jd016542 (2011).Liu, Z., Wang, L. & Wang, S. Comparison of Different GPP Models in China Using MODIS Image and ChinaFLUX Data. Remote Sensing 6, 10215–10231, https://doi.org/10.3390/rs61010215 (2014).Article 
ADS 

Google Scholar 
Edreira, J., Guilpart, N., Sadras, V., Cassman, K. G. & Grassini, P. Water productivity of rainfed maize and wheat: A local to global perspective. Agricultural and Forest Meteorology 259, 364–373, https://doi.org/10.1016/j.agrformet.2018.05.019 (2018).Article 
ADS 

Google Scholar 
Li, H. et al. Water Use Characteristics of Maize-Green Manure Intercropping Under Different Nitrogen Application Levels in the Oasis Irrigation Area Scientia Agricultura Sinica 54, 2608–2618 (2021).
Google Scholar 
Wang, S., Ibrom, A., Bauer-Gottwein, P. & Garcia, M. Incorporating diffuse radiation into a light use efficiency and evapotranspiration model: An 11-year study in a high latitude deciduous forest. Agricultural and Forest Meteorology https://doi.org/10.1016/j.agrformet.2017.10.023 (2018).Article 

Google Scholar 
Cheng, M. High-resolution crop yield and water productivity dataset generated using random forest and remote sensing. Zenodo https://doi.org/10.5281/zenodo.6444614 (2022). More

Maxwell, S. L., Fuller, R. A., Brooks, T. M. & Watson, J. E. M. Biodiversity: The ravages of guns, nets and bulldozers. Nature 536, 143–145 (2016).ADS 
CAS 
PubMed 

Google Scholar 
Uchida, K. & Ushimaru, A. Biodiversity declines due to abandonment and intensification of agricultural lands: Patterns and mechanisms. Ecol. Monogr. 84, 637–658 (2014).
Google Scholar 
Habel, J. C. et al. Butterfly community shifts over two centuries: Shifts in butterfly communities. Conserv. Biol. 30, 754–762 (2016).PubMed 

Google Scholar 
Seibold, S. et al. Arthropod decline in grasslands and forests is associated with landscape-level drivers. Nature 574, 671–674 (2019).ADS 
CAS 
PubMed 

Google Scholar 
Hallmann, C. A. et al. More than 75 percent decline over 27 years in total flying insect biomass in protected areas. PLoS One 12, e0185809 (2017).PubMed 
PubMed Central 

Google Scholar 
Wenzel, M., Schmitt, T., Weitzel, M. & Seitz, A. The severe decline of butterflies on western German calcareous grasslands during the last 30 years: A conservation problem. Biol. Cons. 128, 542–552 (2006).
Google Scholar 
Biesmeijer, J. C. et al. Parallel declines in pollinators and insect-pollinated plants in Britain and the Netherlands. Science 313, 351–354 (2006).ADS 
CAS 
PubMed 

Google Scholar 
Hallmann, C. A., Foppen, R. P. B., van Turnhout, C. A. M., de Kroon, H. & Jongejans, E. Declines in insectivorous birds are associated with high neonicotinoid concentrations. Nature 511, 341–343 (2014).ADS 
CAS 
PubMed 

Google Scholar 
Møller, A. P. Parallel declines in abundance of insects and insectivorous birds in Denmark over 22 years. Ecol. Evol. 9, 6581–6587 (2019).PubMed 
PubMed Central 

Google Scholar 
Wagner, D. L. Insect declines in the anthropocene. Annu. Rev. Entomol. 65, 457–480 (2020).CAS 
PubMed 

Google Scholar 
Habel, J. C., Samways, M. J. & Schmitt, T. Mitigating the precipitous decline of terrestrial European insects: Requirements for a new strategy. Biodivers. Conserv. 28, 1343–1360 (2019).
Google Scholar 
Uhl, B., Wölfling, M. & Fiedler, K. Understanding small-scale insect diversity patterns inside two nature reserves: The role of local and landscape factors. Biodivers. Conserv. 29, 2399–2418 (2020).
Google Scholar 
Stevens, C. J., Dise, N. B., Mountford, J. O. & Gowing, D. J. Impact of nitrogen deposition on the species richness of grasslands. Science 303, 1876–1879 (2004).ADS 
CAS 
PubMed 

Google Scholar 
Thomas, J. A. Butterfly communities under threat. Science 353, 216–218 (2016).ADS 
CAS 
PubMed 

Google Scholar 
Sanders, J. & Hess, J. Benefits of organic farming to environment and society. Thünen Report 65, 362 (2019).
Google Scholar 
Brühl, C. A. & Zaller, J. G. Biodiversity decline as a consequence of an inappropriate environmental risk assessment of pesticides. Front. Environ. Sci. 7, 177 (2019).
Google Scholar 
Brühl, C. A. et al. Direct pesticide exposure of insects in nature conservation areas in Germany. Sci. Rep. 11, 24144 (2021).ADS 
PubMed 
PubMed Central 

Google Scholar 
Wagner, D. L., Grames, E. M., Forister, M. L., Berenbaum, M. R. & Stopak, D. Insect decline in the Anthropocene: Death by a thousand cuts. Proc. Natl. Acad. Sci. USA 118, e2023989118 (2021).CAS 
PubMed 
PubMed Central 

Google Scholar 
Den Boer, P. J. & van Dijk, T. S. Carabid Beetles in A Changing Environment (Agricultural Univ, 1995).
Google Scholar 
Cristescu, M. E. From barcoding single individuals to metabarcoding biological communities: Towards an integrative approach to the study of global biodiversity. Trends Ecol. Evol. 29, 566–571 (2014).PubMed 

Google Scholar 
Hausmann, A. et al. Toward a standardized quantitative and qualitative insect monitoring scheme. Ecol. Evol. 10, 4009–4020 (2020).PubMed 
PubMed Central 

Google Scholar 
Ratnasingham, S. & Hebert, P. D. N. A DNA-based registry for all animal species: The Barcode Index Number (BIN) system. PLoS One 8, e66213 (2013).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Hausmann, A. et al. Genetic patterns in european geometrid moths revealed by the Barcode Index Number (BIN) system. PLoS One 8, e84518 (2013).ADS 
PubMed 
PubMed Central 

Google Scholar 
Padial, J. M., Miralles, A., De la Riva, I. & Vences, M. The integrative future of taxonomy. Front. Zool. 7, 1–14 (2010).
Google Scholar 
Schlick-Steiner, B. C. et al. Integrative taxonomy: A multisource approach to exploring biodiversity. Ann. Rev. Entomol. 55, 421–438 (2010).CAS 

Google Scholar 
Schlick‐Steiner, B. C., Arthofer, W., & Steiner, F. M. Take up the challenge! Opportunities for evolution research from resolving conflict in integrative taxonomy (2014).Fujita, M. K., Leaché, A. D., Burbrink, F. T., McGuire, J. A. & Moritz, C. Coalescent-based species delimitation in an integrative taxonomy. Trends Ecol. Evol. 27, 480–488 (2012).PubMed 

Google Scholar 
Morinière, J. et al. A DNA barcode library for 5,200 German flies and midges (Insecta: Diptera) and its implications for metabarcoding-based biomonitoring. Mol. Ecol. Res. 19, 900–928 (2019).
Google Scholar 
Kortmann, M. et al. Arthropod dark taxa provide new insights into diversity responses to bark beetle infestations. Ecol. Appl. 32, e2516 (2022).PubMed 

Google Scholar 
Porter, T. M. & Hajibabaei, M. Automated high throughput animal CO1 metabarcode classification. Sci. Rep. 8, 1–10 (2018).
Google Scholar 
Boggs, C. L. & Inouye, D. W. A single climate driver has direct and indirect effects on insect population dynamics: Climate drivers of population dynamics. Ecol. Lett. 15, 502–508 (2012).PubMed 

Google Scholar 
Conrad, K. F., Fox, R. & Woiwod, I. P. Monitoring biodiversity: Measuring long-term changes in insect abundance. In Insect Conservation Biology (eds Stewart, A. J. A. et al.) 203–225 (CABI, 2007). https://doi.org/10.1079/9781845932541.0203.Chapter 

Google Scholar 
Flohre, A. et al. Agricultural intensification and biodiversity partitioning in European landscapes comparing plants, carabids, and birds. Ecol. Appl. Publ. Ecol. Soc. Am. 21, 1772–1781 (2011).
Google Scholar 
Emmerson, M. et al. How agricultural intensification affects biodiversity and ecosystem services. In Advances in Ecological Research, vol ***55 43–97 (Elsevier, 2016).
Google Scholar 
Segerer, A. H. & Rosenkranz, E. Das grosse Insektensterben: Was es Bedeutet und was Wir Jetzt tun Müssen (Oekom Verlag, 2019).
Google Scholar 
Batáry, et al. The former Iron Curtain still drives biodiversity-profit trade-offs in German agriculture. Nat. Ecol. Evol. 1, 1279–1284 (2017).PubMed 

Google Scholar 
Kuussaari, M. et al. Extinction debt: A challenge for biodiversity conservation. Trends Ecol. Evol. 24, 564–571 (2009).PubMed 

Google Scholar 
Birkhofer, K., Smith, H. G., Weisser, W. W., Wolters, V. & Gossner, M. M. Land-use effects on the functional distinctness of arthropod communities. Ecography 38, 889–900 (2015).
Google Scholar 
Tscharntke, T., Klein, A. M., Kruess, A., Steffan-Dewenter, I. & Thies, C. Landscape perspectives on agricultural intensification and biodiversity—ecosystem service management. Ecol. Lett. 8, 857–874 (2005).
Google Scholar 
Habel, J. C., Seibold, S., Ulrich, W. & Schmitt, T. Seasonality overrides differences in butterfly species composition between natural and anthropogenic forest habitats. Anim. Conserv. 21, 405–413 (2018).
Google Scholar 
Schmitt, T., Ulrich, W., Delic, A., Teucher, M. & Habel, J. C. Seasonality and landscape characteristics impact species community structure and temporal dynamics of East African butterflies. Sci. Rep. 11, 15103 (2021).CAS 
PubMed 
PubMed Central 

Google Scholar 
Ssymank, A. et al. Praktische Hinweise und Empfehlungen zur Anwendung von Malaisefallen für Insekten in der Biodiversitätserfassung und im Monitoring. Entomol. Verein Krefeld 1, 1–12 (2018).
Google Scholar 
Elbrecht, V., Peinert, B. & Leese, F. Sorting things out: Assessing effects of unequal specimen biomass on DNA metabarcoding. Ecol. Evol. 7, 6918–6926 (2017).PubMed 
PubMed Central 

Google Scholar 
Elbrecht, V. & Steinke, D. Scaling up DNA metabarcoding for freshwater macrozoobenthos monitoring. Freshw. Biol. 64, 380–387 (2019).CAS 

Google Scholar 
Boetzl, F. A. et al. A multitaxa assessment of the effectiveness of agri-environmental schemes for biodiversity management. Proc. Natl. Acad. Sci. 118, 25 (2021).
Google Scholar 
Uhler, J. et al. Relationship of insect biomass and richness with land use along a climate gradient. Nat. Commun. 12, 1–9 (2021).
Google Scholar 
Leray, M. et al. A new versatile primer set targeting a short fragment of the mitochondrial COI region for metabarcoding metazoan diversity: Application for characterizing coral reef fish gut contents. Front. Zool. 10, 34 (2013).PubMed 
PubMed Central 

Google Scholar 
Morinière, J. et al. Species identification in malaise trap samples by DNA barcoding based on NGS Technologies and a scoring matrix. PLoS One 11, e0155497 (2016).PubMed 
PubMed Central 

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

Google Scholar 
Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 17, 10 (2011).
Google Scholar 
Ondov, B. D., Bergman, N. H. & Phillippy, A. M. Interactive metagenomic visualization in a Web browser. BMC Bioinform. 12, 385 (2011).
Google Scholar  More

Goodale, E., Beauchamp, G. & Ruxton, G. D. Mixed-Species Groups of Animals: Behavior, Community Structure, and Conservation (Academic Press, 2017).
Google Scholar 
Krause, J. & Ruxton, G. D. Living in Groups (Oxford University Press, 2002).
Google Scholar 
Stensland, E., Angerbjorn, A. & Berggren, P. Mixed species groups in mammals. Mamm. Rev. 33, 205–223 (2003).
Google Scholar 
Anderson, T. M. et al. Landscape-scale analyses suggest both nutrient and antipredator advantages to Serengeti herbivore hotspots. Ecology 91, 1519–1529 (2010).PubMed 

Google Scholar 
Sinclair, A. R. E. Does interspecific competition or predation shape the African ungulate community? J. Anim. Ecol. 54, 899–918 (1985).
Google Scholar 
Kiffner, C., Kioko, J., Leweri, C. & Krause, S. Seasonal patterns of mixed species groups in large East African mammals. PLoS ONE 9, e113446 (2014).ADS 
PubMed 
PubMed Central 

Google Scholar 
Meise, K., Franks, D. W. & Bro-Jørgensen, J. Using social network analysis of mixed species groups in African savannah herbivores to assess how community structure responds to environmental change. Philos. Trans. R. Soc. B Biol. Sci. 374, 20190009 (2019).
Google Scholar 
de Boer, W. F. & Prins, H. H. T. Large herbivores that thrive mightily but eat and drink as friends. Oecologia 82, 264–274 (1990).ADS 
PubMed 

Google Scholar 
Beaudrot, L., Palmer, M. S., Anderson, T. M. & Packer, C. Mixed-species groups of Serengeti grazers: A test of the stress gradient hypothesis. Ecology. https://doi.org/10.1002/ecy.3163 (2020).Article 
PubMed 

Google Scholar 
He, Q., Bertness, M. D. & Altieri, A. H. Global shifts towards positive species interactions with increasing environmental stress. Ecol. Lett. 16, 695–706 (2013).PubMed 

Google Scholar 
Bertness, M. D. & Callaway, R. Positive interactions in communities. Trends Ecol. Evol. 9, 191–193 (1994).CAS 
PubMed 

Google Scholar 
Fugère, V. et al. Testing the stress-gradient hypothesis with aquatic detritivorous invertebrates: Insights for biodiversity-ecosystem functioning research. J. Anim. Ecol. 81, 1259–1267 (2012).PubMed 

Google Scholar 
Bakker, E. S., Dobrescu, I., Straile, D. & Holmgren, M. Testing the stress gradient hypothesis in herbivore communities: Facilitation peaks at intermediate nutrient levels. Ecology 94, 1776–1784 (2013).PubMed 

Google Scholar 
Hopcraft, J. G. C., Olff, H. & Sinclair, A. R. E. Herbivores, resources and risks: Alternating regulation along primary environmental gradients in savannas. Trends Ecol. Evol. 25, 119–128 (2010).PubMed 

Google Scholar 
Sih, A. Optimal behavior: Can foragers balance two conflicting demands? Science 210, 1041–1043 (1980).ADS 
CAS 
PubMed 

Google Scholar 
Creel, S. & Christianson, D. Relationships between direct predation and risk effects. Trends Ecol. Evol. 23, 194–201 (2008).PubMed 

Google Scholar 
Zollner, P. A. & Lima, S. L. Towards a behavioral ecology of ecological landscapes. Trends Ecol. Evol. 11, 131–135 (1996).PubMed 

Google Scholar 
Brown, J. S., Laundré, J. W. & Gurung, M. The ecology of fear: Optimal foraging, game theory, and trophic interactions. J. Mammal. 80, 385–399 (1999).
Google Scholar 
Gaynor, K. M., Brown, J. S., Middleton, A. D., Power, M. E. & Brashares, J. S. Landscapes of fear: Spatial patterns of risk perception and response. Trends Ecol. Evol. 34, 355–368 (2019).PubMed 

Google Scholar 
Creel, S., Schuette, P. & Christianson, D. Effects of predation risk on group size, vigilance, and foraging behavior in an African ungulate community. Behav. Ecol. 25, 773–784 (2014).
Google Scholar 
Goodale, E., Beauchamp, G., Magrath, R. D., Nieh, J. C. & Ruxton, G. D. Interspecific information transfer influences animal community structure. Trends Ecol. Evol. 25, 354–361 (2010).PubMed 

Google Scholar 
Freeberg, T. M., Eppert, S. K., Sieving, K. E. & Lucas, J. R. Diversity in mixed species groups improves success in a novel feeder test in a wild songbird community. Sci. Rep. 7, 43014 (2017).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Anderson, T. M. et al. The spatial distribution of african savannah herbivores: Species associations and habitat occupancy in a landscape context. Philos. Trans. R. Soc. B Biol. Sci. 371, 20150314 (2016).
Google Scholar 
Arsenault, R. & Owen-Smith, N. Resource partitioning by grass height among grazing ungulates does not follow body size relation. Oikos 117, 1711–1717 (2008).
Google Scholar 
Esmaeili, S. et al. Body size and digestive system shape resource selection by ungulates: A cross-taxa test of the forage maturation hypothesis. Ecol. Lett. 24, 2178–2191 (2021).PubMed 

Google Scholar 
Hopcraft, J. G. C., Anderson, T. M., Pérez-Vila, S., Mayemba, E. & Olff, H. Body size and the division of niche space: Food and predation differentially shape the distribution of Serengeti grazers. J. Anim. Ecol. 81, 201–213 (2012).PubMed 

Google Scholar 
McArthur, C., Banks, P. B., Boonstra, R. & Forbey, J. S. The dilemma of foraging herbivores: Dealing with food and fear. Oecologia 176, 677–689 (2014).ADS 
PubMed 

Google Scholar 
Gagnon, M. & Chew, A. E. Dietary preferences in extant African Bovidae. J. Mammal. 81, 490–511 (2000).
Google Scholar 
Kartzinel, T. R. et al. DNA metabarcoding illuminates dietary niche partitioning by African large herbivores. Proc. Natl. Acad. Sci. U.S.A. 112, 8019–8024 (2015).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Veldhuis, M. P. et al. Cross-boundary human impacts compromise the Serengeti-Mara ecosystem. Science 363, 1424–1428 (2019).ADS 
CAS 
PubMed 

Google Scholar 
Kavwele, C. M. et al. Non-local effects of human activity on the spatial distribution of migratory wildlife in Serengeti National Park, Tanzania. Ecol. Solut. Evid. 3, e12159 (2022).
Google Scholar 
Bijlsma, R. & Loeschcke, V. Environmental stress, adaptation and evolution: An overview. J. Evol. Biol. 18, 744–749 (2005).CAS 
PubMed 

Google Scholar 
Schmitt, M. H., Stears, K. & Shrader, A. M. Zebra reduce predation risk in mixed-species herds by eavesdropping on cues from giraffe. Behav. Ecol. 27, 1073–1077 (2016).
Google Scholar 
Preisser, E. L., Orrock, J. L. & Schmitz, O. J. Predator hunting mode and habitat domain alter nonconsmuptive effects in predator-prey interactions. Ecology 88, 2744–2751 (2007).PubMed 

Google Scholar 
Kiffner, C. et al. Long-term persistence of wildlife populations in a pastoral area. Ecol. Evol. 10, 10000–10016 (2020).PubMed 
PubMed Central 

Google Scholar 
Hopcraft, J. G. C. et al. Competition, predation, and migration: Individual choice patterns of Serengeti migrants captured by hierarchical models. Ecol. Monogr. 84, 355–372 (2014).
Google Scholar 
Fryxell, J. M. Forage quality and aggregation by large herbivores. Am. Nat. 138, 478–498 (1991).
Google Scholar 
Fitzgibbon, C. D. Mixed-species grouping in Thomson’s and Grant’s gazelles: The antipredator benefits. Anim. Behav. 39, 1116–1126 (1990).
Google Scholar 
Brown, J. S. & Kotler, B. P. Hazardous duty pay and the foraging cost of predation. Ecol. Lett. 7, 999–1014 (2004).
Google Scholar 
Stears, K. & Shrader, A. M. Increases in food availability can tempt oribi antelope into taking greater risks at both large and small spatial scales. Anim. Behav. 108, 155–164 (2015).
Google Scholar 
Creel, S. Toward a predictive theory of risk effects: Hypotheses for prey attributes and compensatory mortality. Ecology 92, 2190–2195 (2011).PubMed 

Google Scholar 
Périquet, S. et al. Effects of lions on behaviour and endocrine stress in plains zebras. Ethology 123, 667 (2017).
Google Scholar 
Stears, K., Schmitt, M. H., Wilmers, C. C. & Shrader, A. M. Mixed-species herding levels the landscape of fear. Proc. R. Soc. B Biol. Sci. 287, 20192555 (2020).
Google Scholar 
Schmitt, M. H., Stears, K., Wilmers, C. C. & Shrader, A. M. Determining the relative importance of dilution and detection for zebra foraging in mixed-species herds. Anim. Behav. 96, 151–158 (2014).
Google Scholar 
Meise, K., Franks, D. W. & Bro-Jørgensen, J. Alarm communication networks as a driver of community structure in African savannah herbivores. Ecol. Lett. 23, 293–304 (2020).PubMed 

Google Scholar 
Codron, D., Hofmann, R. R. & Clauss, M. Morphological and physiological adaptations for browsing and grazing. In The Ecology of Browsing and Grazing II (eds Gordon, I. J. & Prins, H. H. T.) 81–125 (Springer, 2019).
Google Scholar 
Odadi, W. O., Karachi, M. K., Abdulrazak, S. A. & Young, T. P. African wild ungulates compete with or facilitate cattle depending on season. Science 333, 1753–1755 (2011).ADS 
CAS 
PubMed 

Google Scholar 
Maestre, F. T., Callaway, R. M., Valladares, F. & Lortie, C. J. Refining the stress-gradient hypothesis for competition and facilitation in plant communities. J. Ecol. 97, 199–205 (2009).
Google Scholar 
de Jonge, M. M. J. et al. Conditional love? Co-occurrence patterns of drought-sensitive species in European grasslands are consistent with the stress-gradient hypothesis. Glob. Ecol. Biogeogr. 30, 1609–1620 (2021).PubMed 
PubMed Central 

Google Scholar 
Franks, D. W., Weiss, M. N., Silk, M. J., Perryman, R. J. Y. & Croft, D. P. Calculating effect sizes in animal social network analysis. Methods Ecol. Evol. 12, 33–41 (2021).
Google Scholar 
Estes, J. A. et al. Trophic downgrading of planet earth. Science 333, 301–306 (2011).ADS 
CAS 
PubMed 

Google Scholar 
Meise, K., Franks, D. W. & Bro-Jørgensen, J. Multiple adaptive and non-adaptive processes determine responsiveness to heterospecific alarm calls in African savannah herbivores. Proc. R. Soc. B Biol. Sci. 285, 20172676 (2018).
Google Scholar 
Blumstein, D. T., Bitton, A. & DaVeiga, J. How does the presence of predators influence the persistence of antipredator behavior? J. Theor. Biol. 239, 460–468 (2006).ADS 
MathSciNet 
PubMed 
MATH 

Google Scholar 
Riggio, J. et al. Lion populations may be declining in Africa but not as Bauer et al. suggest. Proc. Natl. Acad. Sci. 113, 201521506 (2015).
Google Scholar 
Bauer, H. et al. Lion (Panthera leo) populations are declining rapidly across Africa, except in intensively managed areas. Proc. Natl. Acad. Sci. 112, 14894–14899 (2015).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Pettorelli, N., Bro-Jørgensen, J., Durant, S. M., Blackburn, T. & Carbone, C. Energy availability and density estimates in African ungulates. Am. Nat. 173, 698–704 (2009).PubMed 

Google Scholar 
Haile, G. G. et al. Projected impacts of climate change on drought patterns over East Africa. Earth’s Future 8, 1–23 (2020).
Google Scholar 
Devine, A. P., McDonald, R. A., Quaife, T. & Maclean, I. M. D. Determinants of woody encroachment and cover in African savannas. Oecologia 183, 939–951 (2017).ADS 
PubMed 
PubMed Central 

Google Scholar 
Kiffner, C. et al. Long-term population dynamics in a multi-species assemblage of large herbivores in East Africa. Ecosphere 8, e02027 (2017).
Google Scholar 
Prins, H. H. T. & Loth, P. E. Rainfall patterns as background to plant phenology in northern Tanzania. J. Biogeogr. 15, 451–463 (1988).
Google Scholar 
Beattie, K., Olson, E. R., Kissui, B., Kirschbaum, A. & Kiffner, C. Predicting livestock depredation risk by African lions (Panthera leo) in a multi-use area of northern Tanzania. Eur. J. Wildl. Res. 66, 11 (2020).
Google Scholar 
Kasozi, H. & Montgomery, R. A. Variability in the estimation of ungulate group sizes complicates ecological inference. Ecol. Evol. 10, 6881–6889 (2020).PubMed 
PubMed Central 

Google Scholar 
USGS. MOD13Q1 v006 MODIS/Terra Vegetation Indices 16-Day L3 Global 250 m SIN Grid. 10.5067/MODIS/MOD13Q1.006 (2020).R Core Team. R: A Language and Environment for Statistical Computing. http://www.r-project.org/. Accessed January 02, 2022 (2021).Dice, L. R. Measures of the amount of ecologic association between species. Ecology 26, 297–302 (1945).
Google Scholar 
Croft, D. P., James, R. & Krause, J. Exploring Animal Social Networks (Princeton University Press, 2008).
Google Scholar 
Besag, J. & Clifford, P. Generalized Monte Carlo significance tests. Biometrika 76, 633–642 (1989).MathSciNet 
MATH 

Google Scholar 
Hayward, M. W. & Kerley, G. I. H. Prey preferences of the lion (Panthera leo). J. Zool. 267, 309–322 (2005).
Google Scholar 
Codron, D. et al. Diets of savanna ungulates from stable carbon isotope composition of faeces. J. Zool. 273, 21–29 (2007).
Google Scholar 
Kartzinel, T. R. & Pringle, R. M. Multiple dimensions of dietary diversity in large mammalian herbivores. J. Anim. Ecol. 89, 1482–1496 (2020).PubMed 

Google Scholar 
Prins, H. H. T. & Douglas-Hamilton, I. Stability in a multi-species assemblage of large herbivores in East Africa. Oecologia 83, 392–400 (1990).ADS 
CAS 
PubMed 

Google Scholar 
Tournier, E. et al. Differences in diet between six neighbouring groups of vervet monkeys. Ethology 120, 471–482 (2014).
Google Scholar 
Humphries, B. D., Ramesh, T. & Downs, C. T. Diet of black-backed jackals (Canis mesomelas) on farmlands in the KwaZulu-Natal Midlands, South Africa. Mammalia 80, 405–412 (2016).
Google Scholar  More

Thermogravimetric/differential thermogravimetry analyses of corn stoverThermogravimetric/Differential Thermogravimetry (TG/DTG) curves are shown in Fig. 2. The pyrolysis process of corn stover could be divided into three stages. The first stage was the dehydration stage, which occurred at approximately 55–125 °C, and the weight loss was mainly accounted for by water19. The second stage was the pyrolysis stage, which occurred at approximately 200–400 °C and mainly involved the decomposition of cellulose, hemicellulose and a small amount of lignin. This process involved the generation of CO and CO2 and the breaking of carbonaceous polymer bonds20. In addition, a shoulder peak in the range of 265 to 300 °C in the DTG diagram could be caused by side chain decomposition and glycosidic bond cleavage of xylan during the pyrolysis of corn stover21. The third stage was the carbonization stage, which occurred above 400 °C; this stage mainly involved the decomposition of lignin22,23. The carbonization process was relatively slow after 600 °C; this process was called the passive pyrolysis stage24. In general, the TG loss in the pyrolysis process of corn stover was mainly from the moisture in the biomass sample in the first stage. Hemicellulose and cellulose decomposition occurred in the second stage, and lignin decomposition occurred in the third stage25. In this experiment, the minimum pyrolysis temperature for the preparation of biochar was 400 °C. Therefore, the pyrolysis of biochar was relatively complete.Figure 2TG/DTG curves of corn stover.Full size imageCharacterization of biocharYield and specific surface area analysesThe yield and SBET are presented in Table 2. BC, BC-H and BC-OH represent the origin, acid-modified, and base-modified biochar, respectively. The yield of corn stover biochar exhibited a negative correlation with the temperature and decreased from 39.65 to 28.26% when the pyrolysis temperature increased from 400 to 700 °C. This phenomenon could have occurred due to the loss of more volatile substances and the thermal degradation of lignocellulose with increasing temperature, thus reducing the yield of biochar26,27. The SBET of the original biochar showed little difference below 700 °C but increased significantly at 700 °C. Combined with the SEM analysis (Fig. 3), at low temperatures, more ashes on the surface of biochar could block its pores so that the change in SBET was not obvious. At 700 °C, because the ash content significantly reduced and the pyrolysis was more sufficient, the pores of the biochar were more developed, and the SEBT significantly increased. The SBET of the acid/base-modified biochar increased with increasing temperature. The SBET of biochar was larger than that of the original biochar after acid and base modification at 400–600 °C. This phenomenon occurred because the porous structure of biochar was enhanced by acid and base modification28. Moreover, pickling removed most of the inorganic substances in biochar and reduced ash content, while alkali washing removed the tar on the surface of biochar to a certain extent29. However, at 700 °C, the SBET of biochar after acid/base modification was lower than that of the original biochar. Combined with the SEM (Fig. 4), the acid/base modification caused the nanopores of biochar to collapse into mesopores or macropores30. Therefore, the well-developed pore structure of the biochar prepared at 700 °C was destroyed by acid/base modification, resulting in a significant decrease in SBET.Table 2 Yield and SBET of different biochars.Full size tableFigure 3SEM (ZEISS) images of biochar at different pyrolysis temperatures: (a) C1, (b) C8, (c) C12, and (d) C16.Full size imageFigure 4SEM (OPTON) images of C16 biochar and its acid/base modification: (a) C16, (b) C16-H, and (c) C–OH.Full size imageScanning electron microscopy analysisThe C1, C8, C12 and C16 biochars had the highest Cd2+ removal rates at 400, 500, 600 and 700 °C, respectively. Therefore, these BCs were selected for SEM analysis. Figure 3 clearly showed that as the pyrolysis temperature increased from 400 to 700 °C, the pore structure of biochar became more developed, with a smaller pore size and more pores. Although there were numerous pores at 500 °C, the pores were not fully developed and were blocked inside. At 700 °C, the skeleton structure appeared, and the particle size of ash blocked in the pores decreased.By taking C16 biochar with the highest removal rate of Cd2+ as the research object, the changes in the biochar surface before and after modification were compared. C16-H and C16-OH represent acid-modified and base-modified biochar, respectively. After acid/base modification, the ash content on the surface of the biochar decreased, and the pore size increased (Fig. 4). Therefore, some skeleton structures could collapse after corrosion, which was consistent with the previous SBET results. Sun et al. discovered that citric acid-modified biochar would lead to micropore wall collapse and micropore loss, resulting in a reduction in SBET31. This finding was in agreement with the results of our study.Fourier transform infrared spectroscopy analysisThe FTIR spectra of biochar at different pyrolysis temperatures are presented in Fig. 5a.Figure 5FTIR spectra of corn stover biochar: (a) different pyrolysis temperatures and (b) different modification treatments.Full size imageAs the pyrolysis temperature increased from 300 to 700 °C, the absorption peak intensity showed a downwards trend. There was a remarkable decrease in features associated with stretch O–H (3400 cm−1)32. The vibration peaks of C–H (2924 cm−1) and C=O (1610 cm−1) decreased with increasing temperature, which could be due to the reduction in –CH2 and –CH3 groups of small molecules and the pyrolysis of C=O into gas or liquid byproducts at high temperatures33. In addition, the peak at 1435 cm−1 was identified as the vibration of C=C bonds belonging to the aromatic skeleton of biochar. A decrease in the absorbance peaks was found at 1115 cm−1, which corresponded to C–O–C bonds. The ratio of intensities for C=C/C=O (1550–1650 cm−1) and C–O–C (1115 cm−1) to the shoulder (1100–1200 cm−1) gradually decreased, and the loss of –OH at 3444 cm−1 indicated that the oxygen content in biochar reduced. The cellulose and wood components were dehydrated, and the degree of biochar condensation increased at higher temperatures. The bending vibration peaks of Ar–H at 856 and 877 cm−1 changed little at different temperatures, which showed that the aromatic rings were relatively stable below 700 °C34. Combined with the above analysis the condensation degree of biochar increased gradually above 400 °C35,36. In summary, as the pyrolysis temperature increased, the degree of aromatization of biochar improved, and the numbers of oxygen-containing functional groups decreased continuously.Figure 5b showed that after acid/base modification, the absorbance peaks at 3444 cm−1, 1610 cm−1 and 1115 cm−1 increased, indicating that the number of oxygen-containing functional groups increased. However, the stretching vibration peak of aromatic ring skeleton C=C (1435 cm−1) and the bending vibration peaks of Ar–H (856–877 cm−1) changed little. The number of functional groups of acid-modified biochar increased more than that of alkali-modified biochar. Mahdi et al. found that acid modification increased the number of functional groups in a study of biochar modification37. After acid/base modification, the number of oxygen-containing functional groups, such as hydroxyl and carboxyl groups, increased.Optimization of biocharFigure 6 illustrates that the removal rates of Cd2+ by corn stover biochar (original, acid-modified, and base-modified biochars) consistently increased with increasing pyrolysis temperature. The highest removal rate reached 95.79% at 700 °C. The removal rate decreased after modification, especially after pickling. The results showed that C16 biochar had the best removal effect on Cd2+.Figure 6Cd2+ removal rate of different biochars (BC: original biochar, BC-OH: alkali-modified biochar, and BC-H: acid-modified biochar).Full size imageIntuitive and variance analyses were employed to explore the influences of biochar preparation conditions on the removal rate of Cd2+.

1.

Intuitive analysis
The intuitive analysis of the orthogonal experiment is shown in Table 3 and Fig. 7. The pyrolysis temperature had the most significant influence on the removal of Cd2+, followed by the retention time and finally the heating rate. Therefore, the optimal conditions for biochar preparation were a pyrolysis temperature of 700 °C, a retention time of 2.5 h, and a heating rate of 5 °C/min.

2.

Variance analysis
Variance analysis showed that the effect of pyrolysis temperature on the removal rate of Cd2+ was very significant (Table 4). The effects of retention time and heating rate were not significant. This phenomenon was consistent with the conclusions obtained in the intuitive analysis.

Table 3 Intuitive analyses of influencing factors of biochar preparation.Full size tableFigure 7Intuitive analysis diagram of influencing factors for biochar preparation.Full size imageTable 4 Variance analysis.Full size tableAnalysis of adsorption mechanismThe SBET of the unmodified biochar did not change significantly with temperature, which indicated that SBET could potentially not be a critical factor for Cd2+ adsorption. Qi et al. obtained a similar conclusion when studying the adsorption of Cd2+ in water by chicken litter biochar38. In addition to SBET, the four primary mechanisms involved in the removal of heavy metal ions by biochar were as follows: (1) Ion exchange: the alkali or alkaline earth metals in biochar (K+, Ca2 +, Na+, and Mg2+) were the dominant cations in ion exchange39. (2) The complexation of oxygen-containing functional groups mainly included hydroxyl and carboxyl groups40. (3) Mineral precipitation: Cd2+ was precipitated by minerals on the surface of biochar to form Cd3(PO4)2 and CdCO341. Soluble cadmium precipitated with some anions released by biochar, such as CO32−, PO43− and OH−42,43. (4) π electron interaction: Cd2+ coordinated with the π electrons of C=C or C=O at low pyrolysis temperatures43,44. Biochar contains more aromatic structures at high pyrolysis temperatures, which could provide more π electrons. Therefore, the π electron interaction in adsorption of Cd2+ was effectively enhanced45.C1, C8, C12 and C16 were selected to study the adsorption mechanism. Related physicochemical properties are given in Table 5.Table 5 Physicochemical properties of biochar at different pyrolysis temperatures.Full size tableThe CEC of biochar gradually increased as the pyrolysis temperature increased, reaching a maximum at 600 °C and slightly decreasing at 700 °C. This phenomenon could have occurred because the crystalline minerals under high pyrolysis temperatures inhibited the exchange of cations on the surface of biochar with Cd2+ in aqueous solution46. Nevertheless, CEC did not change significantly with temperature; thus, CEC was not the main adsorption mechanism. With increasing pyrolysis temperature, the number of acidic functional groups decreased gradually, while the number of alkaline functional groups increased. The main functional groups used to remove Cd2+ were generally considered acidic oxygen-containing functional groups. However, the number of these functional groups decreased with increasing pyrolysis temperature, which weakened the complexation on the surface of the biochar. However, this result was contradictory to the results of Cd2+ adsorption. Therefore, the functional groups were not the main adsorption mechanism.To further explore the adsorption mechanism of Cd2+, the biochar before and after the adsorption of Cd2+ was characterized by XRD. As shown in Fig. 7a, C16-100Cd and C16-200Cd represented the biochar after Cd2+ adsorption when the concentrations of cadmium solution were 100 mg/l and 200 mg/l, respectively. The results showed that new peaks appeared at 30.275° and 36.546° after adsorption, corresponding to CdCO3. The spike at 29.454° was due to Cd(OH)2. Additionally, the intensity of the CdCO3 peak increased significantly from C16-100Cd to C16-200Cd, indicating that mineral precipitation occurred in adsorption. Liu et al. found similar results in a study on removing Cd2+ from water by blue algae biochar12. However, as the concentration of Cd2+ increased from 0 to 200 mg/L, the diffraction peak at 2θ = 29.454° first increased and then decreased. This because the peak position of CaCO3 at 2θ = 29.369° was very close to Cd(OH)2 at 2θ = 29.454°. At low concentrations, the production of Cd(OH)2 was greater than that of CdCO3. When the initial concentration of Cd2+ increased, more CO32− released by CaCO3 combined with Cd2+ to form CdCO3, resulting in a reduction in the diffraction peak.As presented in Fig. 8b, the peak intensities of CdCO3 and Cd(OH)2 gradually increase with increasing pyrolysis temperature. On the one hand, this phenomenon could be ascribed to the increase in the mineral content of biochar with increasing pyrolysis temperature. On the other hand, the pH value of biochar increased with increasing pyrolysis temperature. In this way, more OH− was released, thus forming more Cd(OH)2. Wang et al. obtained similar results42. Moreover, the peak intensity of KCl at 2θ = 28.347° decreased after adsorption, as shown in Fig. 8a, which indicated that ion exchange took part in adsorption.Figure 8XRD images: (a) before and after adsorption of Cd2+ on C16 biochar and (b) Cd2+ adsorption by biochar at different pyrolysis temperatures.Full size imageIn addition, the FTIR spectra showed that the number of functional groups, such as C=C and C=O, in biochar decreased with increasing pyrolysis temperature, leading to the weakening of cation–π interactions between Cd2+ and C=C and C=O. In contrast, due to the enhanced aromatization of functional groups on the surface of biochar, many lone pair electrons existed in the electron-rich domains of the graphene-like structure, which in turn enhanced the cation–π interactions. Harvey et al., based on the study of Cd2+ adsorption by plant biochar, concluded that the electron-rich domain bonding mechanism between Cd2+ and the graphene-like structure on the surface of biochar played a more significant role in biochar with a high degree of carbonization45. Therefore, π-electron interactions could play a dominant role in Cd2+ adsorption on high-temperature pyrolysis biochar. Moreover, the results showed that the number of alkaline functional groups increased while acidic functional groups decreased with the increase in pyrolysis temperature. It is generally believed that acidic functional groups could withdraw electrons, and basic functional groups could donate electrons47,48. The biochar with higher pyrolysis temperature contained more alkaline functional groups, which improved the electron donating ability of biochar and enhanced the cation–π electron effect.In summary, mineral precipitation and π electron coordination were the main mechanisms of removing Cd2+ from water by corn stover biochar. This phenomenon explained why the Cd2+ removal rate of acid/base–modified biochar decreased. After modification, the functional groups on the surface of biochar increased, but the inorganic minerals were removed. Pickling resulted in the loss of soluble minerals and alkaline functional groups on the surface of biochar, which was not conducive to adsorption49. After alkaline washing, more PO43−, CO32− and HCO3− were released, thereby reducing the mineral precipitation50,51. Since NaOH had a weaker destructive effect than HCl and introduced some OH−, alkaline washing had little effect on the removal rate of Cd2+.Adsorption isotherm and adsorption kineticsAdsorption isothermThe adsorption isotherms were fitted with Langmuir (Eq. 3) and Freundlich (Eq. 4) models, as shown in Fig. 9, and the fitting parameters are listed in Table 6.Figure 9Adsorption isotherm.Full size imageTable 6 Fitting parameters of the adsorption isotherm model.Full size tableThe Langmuir model (R2  > 0.963) was more suitable than the Freundlich model (R2  > 0.919), indicating that the adsorption sites of biochar were evenly distributed, and adsorption was mainly monolayer. Parameter KL reflected the difficulty of adsorption and was generally divided into four types: unfavourable (KL  > 1), favourable (0  More