Philippa Kaur

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

Bates, H. W. XXXII. Contributions to an insect fauna of the Amazon Valley. Lepidoptera: Heliconidæ. Trans. Linn. Soc. Lond 23, 495–566 (1862).Article 

Google Scholar 
Müller, F. Ituna and thyridia: A remarkable case of mimicry in butterflies. Trans. Entomol. Soc. Lond. 1879, 20–29 (1879).
Google Scholar 
Baxter, S. W. et al. Convergent evolution in the genetic basis of Müllerian mimicry in Heliconius butterflies. Genetics 180, 1567–1577 (2008).Article 
CAS 

Google Scholar 
Sheppard, P. M., Turner, J. R. G., Brown, K., Benson, W. & Singer, M. Genetics and the evolution of Muellerian mimicry in Heliconius butterflies. Philos. Trans R. Soc. Lond. B, Biol. Sci. 308, 433–610 (1985).Article 
ADS 

Google Scholar 
Mallet, J. & Gilbert, L. E. Jr. Why are there so many mimicry rings? Correlations between habitat, behaviour and mimicry in Heliconius butterflies. Biol. J. Lin. Soc. 55, 159–180 (1995).
Google Scholar 
Brower, A. V. Parallel race formation and the evolution of mimicry in Heliconius butterflies: A phylogenetic hypothesis from mitochondrial DNA sequences. Evolution 50, 195–221 (1996).Article 
CAS 

Google Scholar 
Wilson, J. S. et al. North American velvet ants form one of the world’s largest known Müllerian mimicry complexes. Curr. Biol. 25, R704–R706. https://doi.org/10.1016/j.cub.2015.06.053 (2015).Article 
CAS 
PubMed 

Google Scholar 
Wilson, J. S., Williams, K. A., Forister, M. L., Von Dohlen, C. D. & Pitts, J. P. Repeated evolution in overlapping mimicry rings among North American velvet ants. Nat. Commun. 3, 1272 (2012).Article 
ADS 

Google Scholar 
Wilson, J. S., Pan, A. D., Limb, E. S. & Williams, K. A. Comparison of African and North American velvet ant mimicry complexes: Another example of Africa as the ‘odd man out’. PLoS ONE 13, e0189482. https://doi.org/10.1371/journal.pone.0189482 (2018).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Plowright, R. & Owen, R. E. The evolutionary significance of bumble bee color patterns: A mimetic interpretation. Evolution 34, 622–637 (1980).Article 
CAS 

Google Scholar 
Williams, P. The distribution of bumblebee colour patterns worldwide: Possible significance for thermoregulation, crypsis, and warning mimicry. Biol. J. Lin. Soc. 92, 97–118 (2007).Article 

Google Scholar 
Hines, H. M. & Williams, P. H. Mimetic colour pattern evolution in the highly polymorphic Bombus trifasciatus (Hymenoptera: Apidae) species complex and its comimics. Zool. J. Linn. Soc. 166, 805–826 (2012).Article 

Google Scholar 
Koch, J. B., Rodriguez, J., Pitts, J. P. & Strange, J. P. Phylogeny and population genetic analyses reveals cryptic speciation in the Bombus fervidus species complex (Hymenoptera: Apidae). PLoS ONE 13, e0207080 (2018).Article 

Google Scholar 
Ezray, B. D., Wham, D. C., Hill, C. E. & Hines, H. M. Unsupervised machine learning reveals mimicry complexes in bumblebees occur along a perceptual continuum. Proc. R. Soc. B 286, 20191501 (2019).Article 

Google Scholar 
Bateson, W. The alleged “Aggressive Mimicry” of volucellæ. Nature 46, 585 (1892).Article 
ADS 

Google Scholar 
Poulton, E. B. The volucellœ as alleged examples of variation “almost unique among animals”. Nature 47, 126 (1892).Article 
ADS 

Google Scholar 
Cockerell, T. D. New social bees. Psyche A J. Entomol. 24, 120–128 (1917).Article 

Google Scholar 
Koch, J., Strange, J. & Williams, P. In: Bumble bees of the western United States (US Forest Service, San Francisco California, 2012).
Google Scholar 
Williams, P. H., Thorp, R. W., Richardson, L. L. & Colla, S. R. In: Bumble bees of North America: An identification guide Vol. 87 (Princeton University Press, Princeton, 2014).
Google Scholar 
Ruxton, G. D., Franks, D. W., Balogh, A. C. & Leimar, O. Evolutionary implications of the form of predator generalization for aposematic signals and mimicry in prey. Evol Int. J. Org. Evol. 62, 2913–2921 (2008).Article 

Google Scholar 
Rowe, C., Lindström, L. & Lyytinen, A. The importance of pattern similarity between Müllerian mimics in predator avoidance learning. Proc. R. Soc. Lond. Ser. B Biol. Sci. 271, 407–413 (2004).Article 

Google Scholar 
Beatty, C. D., Beirinckx, K. & Sherratt, T. N. The evolution of Müllerian mimicry in multispecies communities. Nature 431, 63 (2004).Article 
ADS 
CAS 

Google Scholar 
Chittka, L. & Osorio, D. Cognitive dimensions of predator responses to imperfect mimicry. PLoS Biol. 5, e339 (2007).Article 

Google Scholar 
Dittrigh, W., Gilbert, F., Green, P., McGregor, P. & Grewcock, D. Imperfect mimicry: A pigeon’s perspective. Proc. R. Soc. Lond. Ser. B Biol. Sci. 251, 195–200 (1993).Article 
ADS 

Google Scholar 
Sherratt, T. N., Whissell, E., Webster, R. & Kikuchi, D. W. Hierarchical overshadowing of stimuli and its role in mimicry evolution. Anim. Behav. 108, 73–79 (2015).Article 

Google Scholar 
Beatty, C. D., Bain, R. S. & Sherratt, T. N. The evolution of aggregation in profitable and unprofitable prey. Anim. Behav. 70, 199–208 (2005).Article 

Google Scholar 
Kazemi, B., Gamberale-Stille, G., Tullberg, B. S. & Leimar, O. Stimulus salience as an explanation for imperfect mimicry. Curr. Biol. 24, 965–969 (2014).Article 
CAS 

Google Scholar 
Kikuchi, D. W., Dornhaus, A., Gopeechund, V. & Sherratt, T. N. Signal categorization by foraging animals depends on ecological diversity. Elife. 8, e43965 (2019).Article 

Google Scholar 
Rapti, Z., Duennes, M. A. & Cameron, S. A. Defining the colour pattern phenotype in bumble bees (Bombus): A new model for evo devo. Biol. J. Lin. Soc. 113, 384–404 (2014).Article 

Google Scholar 
Wilson, J. S., Sidwell, J. S., Forister, M. L., Williams, K. A. & Pitts, J. P. Thistledown velvet ants in the desert mimicry ring and the evolution of white coloration: Müllerian mimicry, camouflage and thermal ecology. Biol. Lett. 16, 20200242 (2020).Article 

Google Scholar 
Ascher, J. & Pickering, J. Discover Life bee species guide and world checklist (Hymenoptera: Apoidea: Anthophila) (2019).iNaturalist. Available from https://www.inaturalist.org. Accessed [2022].Bombus Latreille, 1802 in GBIF Secretariat (2021). GBIF Backbone Taxonomy. Checklist dataset https://doi.org/10.15468/39omei accessed via GBIF.org on 2021-12-03. More

Longline fishing boats such as these at Honolulu’s harbour in Hawaii must respect a large no-fishing zone off the western side of the archipelago.Credit: Sarah Medoff

Large no-fishing areas can drive the recovery of commercially valuable fish species, a study suggests. Ten years’ worth of fisheries data have shown that catch rates of two important types of tuna increased drastically in the vicinity of a marine protected area surrounding the northwestern Hawaiian islands.“It’s a win–win for fish and fishermen,” says Jennifer Raynor, an economist at the University of Wisconsin–Madison and a co-author of the study, which was published on 20 October in Science1.The results highlight the value of large-scale marine protected areas — a type of environmental management that has emerged in the past two decades, mostly in the Pacific Ocean, says Kekuewa Kikiloi, who studies Hawaiian culture at the University of Hawaii at Mānoa. Countries around the world have committed to protecting 30% of their land and oceans by 2030.Previous research showed that marine protected areas can help to restore populations of creatures that don’t move around much or at all, such as corals2 and lobsters3. Raynor and her colleagues wanted to test whether the areas could also drive the recovery of migratory species and provide spillover benefits for fisheries. The researchers looked at one of the largest such areas in the world, the 1.5-million-square-kilometre Papahānaumokuākea Marine National Monument, which was created in 2006 and expanded in 2016 to protect biological and cultural resources.The team focused on the Hawaiian ‘deep-set’ longline fishery, which mainly targets yellowfin tuna (Thunnus albacares) and bigeye tuna (Thunnus obesus).The researchers analysed catch data collected on fishing vessels between 2010 and late 2019. Then, they compared catch rates at various distances up to 600 nautical miles (1,111 kilometres) from the protected area, before and after its expansion in 2016. (The protected area itself currently extends for 200 nautical miles from the northwestern part of the Hawaiian archipelago.) They found that after the expansion, catch rates — defined as the number of fish caught for every 1,000 hooks deployed — went up, and that the increases were greater the closer the boats were to the no-fishing zone. At distances of up to 100 nautical miles, the catch rate for yellowfin tuna increased by 54%, and that for bigeye tuna by 12%. Some other types of catch rate also increased, but not by equally significant margins.The size of the Papahānaumokuākea Marine National Monument — more than three times the surface area of California — probably played a part in the positive effects, as did its shape. It spans about 2,000 kilometres from west to east, protecting large swathes of ocean waters at tropical latitudes. This means that tropical fish such as yellowfin and bigeye tuna — which tend to move along an east–west axis to stay in their preferred temperature range — can travel a long way and still stay in the no-fishing zone.What’s more, says Raynor, Papahānaumokuākea is a spawning ground for yellowfin tuna. Because the animals don’t travel far from their birthplace, the no-take zone provides refuge from fishing, helping tuna to aggregate and reproduce.“It is exciting to see that there are benefits to the fishing industry from this marine protected area,” says David Kroodsma, director of research and innovation at Global Fishing Watch in Oakland, California, a US non-governmental organization that monitors fishing activity worldwide. However, he adds, it’s unclear whether the results can be generalized to other areas of the world.Regardless, the findings could help others to design marine protected areas so that benefits trickle down to fisheries, says Steve Gaines, a marine ecologist at the University of California, Santa Barbara. The study, he says, “provides a platform to definitively evaluate what is working and what isn’t”.Co-managed by Indigenous populations, the state of Hawaii and the US government, Papahānaumokuākea is an example of a collaborative management strategy that bridges Indigenous knowledge and modern science, Kikiloi says. The approach, he adds, “can work successfully in other places too, if given a chance”. More

We hypothesized that functionally distinct forest types can be mapped at moderate spatial resolutions, using a combination of canopy foliar traits and canopy structure information. Our analysis of LiDAR and imaging spectroscopy data at spatial resolutions ranging from 4 to 200 m (16 m2–40,000 m2), with an emphasis on the 30 m (900 m2) spaceborne hyperspectral spatial resolution, reveals that few remotely sensed canopy properties are needed to successfully identify ecologically distinct forest types at two diverse tropical forest sites in Malaysian Borneo. In testing our second hypothesis that mapped forest types exhibit distinct ecosystem function, we found that forest types identified using remotely sensed leaf P, LMA, Max H, and canopy cover at 20 m height (Cover20) closely align with forest types defined from field-based floristic surveys29,30,31,32,33 and inventory plot-based measurements of growth and mortality rates (Fig. 4b). Our approach, however, enables mapping of their entire spatial extent (Fig. 1) and reveals important structural and functional variation within areas characterized as a single forest type in previous studies (Fig. 3). Current and forthcoming satellite hyperspectral platforms, including PRISMA (30 m), CHIME (20–30 m), and SBG (30 m), have or will have comparable spectral resolution, higher temporal revisits, and much greater geographic coverage. The ability to conduct this type of analysis using remote sensing measurements at 30 m resolution suggests that our method can be applied to these emerging spaceborne imaging spectroscopy data to reveal important differences in structure and function across the world’s tropical forests.Nested functional forest types revealedTo test our first hypothesis, rather than making an a priori decision about the number of k-means clusters (k), we explored the capacity of remotely sensed data to reveal ecologically relevant variation in forest types. Baldeck and Asner took a similar unsupervised approach to estimating beta diversity in South Africa34. Because the choice of k directly influences analysis outcomes, careful selection of k is required. Different approaches for identifying the number of clusters, using the Gapk and Wk elbow metrics35, yielded varying optimal numbers of clusters for the Sepilok and Danum landscapes (Fig. 1, Supplementary Figs. 4 and 5). However, at both sites, a comparison of results based on different values of k revealed ecologically meaningful structural and functional differences and graduated transitions between forest types (Fig. 2, Supplementary Figs. 7 and 8), indicating that the exploration of traits that aggregate or separate forest types as k changes is a valuable exercise. Overlap between the remotely sensed forest type boundaries and inventory plots within distinct forest types indicate that the series of clustered forests align closely with forest types defined based on in situ data on species composition and ecosystem structure. In part, this type of analysis requires careful selection of the number of clusters. Additionally, however, we gained valuable insights via the exploration of varying numbers of clusters as it relates to biologically meaningful categorization of forest types. Extending this method to other parts of the tropics will require similar decision-making, which will either require user input, or the development of robust automated algorithms for selecting k.Forest types capture differences in ecosystem dynamicsWe further evaluated the canopy traits and structural attributes that were most critical for mapping distinct forest types, hypothesizing that mapped forest types exhibit distinct ecosystem function. Forest types revealed by the cluster analyses were distributed along the leaf economic spectrum, where the leaf economic spectrum characterizes a tradeoff in plant growth strategies36. LMA, which can covary strongly with leaf N and P, is a key indicator of plant growth strategies along the spectrum37. At the slow-return end of the leaf economics spectrum, plants in nutrient-poor conditions with low leaf nutrient concentrations invest in leaf structure and defense, expressed as high LMA, strategizing longer-lived, tougher leaves with slower decomposition rates. This strategy comes at the cost of slower growth. At the quick-return end of the spectrum, plants in nutrient-rich environments with higher leaf nutrient concentrations invest less in structure and defense, enabling faster growth and more rapid leaf turnover, i.e., shorter leaf lifespans. This quick-return growth strategy supports higher photosynthetic rates and more rapid carbon gain36.In this study, the principal components and clustering results yielded forest types that are indicative of community level differences associated with leaf economic spectrum differences. The nutrient rich sites (Danum1 and Danum2, Supplementary Fig. 8) show high canopy N and P and low LMA compared to the nutrient poor and acidic sites (Sandstone and Kerangas), which contributes to lower leaf photosynthetic capacity (Vcmax) and growth (Fig. 4b). Foliar N:P also increased with site fertility, confirming that tropical forests are primarily limited by phosphorus, and not nitrogen38,39, with large implications for carbon sequestration in these forests. Orthogonal differences in canopy structure and architecture between Danum forest types and Sepilok Sandstone and Alluvial forests could be indicative of ecosystem scale differences in the sensitivity of these forests to endogenous disturbance processes40.The significant differences in aboveground carbon stocks and growth and mortality rates between forest types further suggests strong differences in ecosystem dynamics. In general, growth rates varied inversely to aboveground carbon, and higher aboveground carbon corresponded to lower mortality rates. As an example, the Sepilok sandstone forests, which are largely comprised of slow-growing dipterocarp species29,33, had the highest median aboveground carbon (236 Mg C ha−1), with higher canopy P and N, and lower LMA. The taller canopy and low canopy leaf nutrient concentrations are consistent with the low growth and mortality rates found in the sandstone forest, indicating a slow-growth strategy yielding larger trees and higher aboveground carbon stocks. In contrast, alluvial forests exhibit high turnover with mortality and growth rates higher relative to Sandstone forests corresponding to lower aboveground carbon on average. Kerangas forests exhibited low aboveground carbon despite an intermediate plot-level growth rate, and mortality rates that were significantly lower than the Danum or alluvial forest types. Kerangas forests, which were characterized by the highest LMA, lowest foliar P and N (Fig. 2a), and the lowest plot-level aboveground carbon density (186 Mg C ha−1; Fig. 4a), are known to have higher stem densities, lower canopy heights, and long-lived leaves5,32,41, suggesting well-developed strategies for nutrient retention42. Interestingly, despite significantly different aboveground carbon and demography, the kerangas and sandstone forests did not differ in LAI or canopy architecture (P:H); although maximum height, Cover20, and Hpeak LAI were significantly higher in the sandstone forest, highlighting the need to account for differences beyond LAI when scaling processes from leaves to ecosystems.In addition, when three forest types were distinguished at Sepilok, the alluvial inventory plot had significantly higher aboveground carbon than the remote sensing-derived alluvial forest extent (Fig. 4a, p  More

Sakschewski, B. et al. Leaf and stem economics spectra drive diversity of functional plant traits in a dynamic global vegetation model. Glob. Change Biol. 21, 2711–2725 (2015).Article 

Google Scholar 
Berzaghi, F. et al. Towards a new generation of trait-flexible vegetation models. Trends Ecol. Evol. 35, 191–205 (2020).Article 
PubMed 

Google Scholar 
Bruelheide, H. et al. Global trait–environment relationships of plant communities. Nat. Ecol. Evol. 2, 1906–1917 (2018).Article 
PubMed 

Google Scholar 
Joswig, J. S. et al. Climatic and soil factors explain the two-dimensional spectrum of global plant trait variation. Nat. Ecol. Evol. 6, 36–50 (2021).Article 
PubMed 
PubMed Central 

Google Scholar 
van Bodegom, P. M., Douma, J. C. & Verheijen, L. M. A fully traits-based approach to modeling global vegetation distribution. Proc. Natl Acad. Sci. USA 111, 13733–13738 (2014).PubMed Central 

Google Scholar 
Moreno Martínez, A. et al. A methodology to derive global maps of leaf traits using remote sensing and climate data. Remote Sens. Environ. 218, 69–88 (2018).Article 

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

Google Scholar 
Kattge, J. et al. TRY—a global database of plant traits. Glob. Change Biol. 17, 2905–2935 (2011).Article 

Google Scholar 
Kattge, J. et al. TRY plant trait database-enhanced coverage and open access. Glob. Change Biol. 26, 119–188 (2020).Article 

Google Scholar 
Jetz, W. et al. Monitoring plant functional diversity from space. Nat. Plants 2, 16024 (2016).Article 
PubMed 

Google Scholar 
Butler, E. E. et al. Mapping local and global variability in plant trait distributions. Proc. Natl Acad. Sci. USA 114, E10937–E10946 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
Boonman, C. C. et al. Assessing the reliability of predicted plant trait distributions at the global scale. Glob. Ecol. Biogeogr. 29, 1034–1051 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
Madani, N. et al. Future global productivity will be affected by plant trait response to climate. Sci. Rep. 8, 2870 (2018).Vallicrosa, H. et al. Global distribution and drivers of forest biome foliar nitrogen to phosphorus ratios (N:P). Glob. Ecol. Biogeogr. 31, 861–871 (2022).Article 

Google Scholar 
Meyer, H. & Pebesma, E. Predicting into unknown space? Estimating the area of applicability of spatial prediction models. Methods Ecol. Evol. 12, 1620–1633 (2021).Article 

Google Scholar 
Schiller, C. et al. Deep learning and citizen science enable automated plant trait predictions from photographs. Sci. Rep. 11, 16395 (2021).Aguirre-Gutiérrez, J. et al. Pantropical modelling of canopy functional traits using sentinel-2 remote sensing data. Remote Sens. Environ. 252, 112–122 (2021).Article 

Google Scholar 
Homolova, L. et al. Review of optical-based remote sensing for plant trait mapping. Ecol. Complex. 15, 1–16 (2013).Article 

Google Scholar 
Van Cleemput, E. et al. The functional characterization of grass-and-shrubland ecosystems using hyperspectral remote sensing: trends, accuracy and moderating variables. Remote Sens. Environ. 209, 747–763 (2018).Article 

Google Scholar 
Kattenborn, T., Fassnacht, F. E. & Schmidtlein, S. Differentiating plant functional types using reflectance: which traits make the difference? Remote Sens. Ecol. Conserv. 5, 5–19 (2019).Article 

Google Scholar 
Hauser, L. T. et al. Explaining discrepancies between spectral and in-situ plant diversity in multispectral satellite earth observation. Remote Sens. Environ. 265, 112684 (2021).Article 

Google Scholar 
Wäldchen, J. & Mäder, P. Plant species identification using computer vision techniques: a systematic literature review. Arch. Comput. Methods Eng. 25, 507–543 (2018).Article 
PubMed 

Google Scholar 
Jones, H. G. What plant is that? Tests of automated image recognition apps for plant identification on plants from the British flora. AoB Plants 12, plaa052 (2020).Article 
PubMed 
PubMed Central 

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

Google Scholar 
WÜest, R. O. et al. Macroecology in the age of big data—where to go from here? J. Biogeogr. 47, 1–12 (2020).Article 

Google Scholar 
Mäder, P. et al. The Flora Incognita app—interactive plant species identification. Methods Ecol. Evol. 12, 1335–1342 (2021).Article 

Google Scholar 
Di Cecco, G. J. et al. Observing the observers: how participants contribute data to iNaturalist and implications for biodiversity science. BioScience 71, 1179–1188 (2021).Article 

Google Scholar 
Mahecha, M. D. et al. Crowd-sourced plant occurrence data provide a reliable description of macroecological gradients. Ecography 44, 1131–1142 (2021).Article 

Google Scholar 
Botella, C. et al. Jointly estimating spatial sampling effort and habitat suitability for multiple species from opportunistic presence-only data. Methods Ecol. Evol. 12, 933–945 (2021).Article 

Google Scholar 
iNaturalist Research-Grade Observations (GBIF, accessed 5 January 2022); https://www.gbif.org/dataset/50c9509d-22c7-4a22-a47d-8c48425ef4a7Callaghan, C. T. et al. Three frontiers for the future of biodiversity research using citizen science data. BioScience 71, 55–63 (2020).
Google Scholar 
Dickinson, J. L., Zuckerberg, B. & Bonter, D. N. Citizen science as an ecological research tool: challenges and benefits. Ann. Rev. Ecol. Evol. Syst. 41, 149–172 (2010).Article 

Google Scholar 
Kosmala, M. et al. Assessing data quality in citizen science. Front. Ecol. Environ. 14, 551–560 (2016).Article 

Google Scholar 
Boakes, E. H. et al. Patterns of contribution to citizen science biodiversity projects increase understanding of volunteers’ recording behaviour. Sci. Rep. 6, 33051 (2016).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Bowler, D.E. et al. Temporal trends in the spatial bias of species occurrence records. Ecography 2022, e06219 (2022). https://doi.org/10.1111/ecog.06219GBIF Occurrence Download (GBIF, 4 January 2022); https://doi.org/10.15468/dl.34tjreBruelheide, H. et al. sPlot—a new tool for global vegetation analyses. journal of vegetation science. J. Veg. Sci. 30, 161–186 (2019).Article 

Google Scholar 
Sabatini, F. et al. sPlotOpen—an environmentally balanced, open access, global dataset of vegetation plots. Glob. Ecol. Biogeogr. 30, 1740–1764 (2021).Article 

Google Scholar 
Whittaker, R.H. et al. Communities and Ecosystems (Macmillan/Collier Macmillan, 1970).Olson, D. M. et al. Terrestrial ecoregions of the world: a new map of life on Earth: a new global map of terrestrial ecoregions provides an innovative tool for conserving biodiversity. BioScience 51, 933–938 (2001).Article 

Google Scholar 
Joswig, J., Wirth, C. & Schuman, M. Climatic and soil factors explain the two-dimensional spectrum of global plant trait variation. Nat. Ecol. Evol. 6, 36–50 (2022).Article 
PubMed 

Google Scholar 
Díaz, S. et al. The global spectrum of plant form and function. Nature 529, 167–171 (2016).Article 
PubMed 

Google Scholar 
Ploton, P. et al. Spatial validation reveals poor predictive performance of large-scale ecological mapping models. Nat. Commun. 11, 4540 (2020).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Meyer, H. & Pebesma, E. Spatial validation reveals poor predictive performance of large-scale ecological mapping models. Methods Ecol. Evol. 12, 1620–1633 (2021).Article 

Google Scholar 
Schrodt, F. et al. Bhpmf—a hierarchical Bayesian approach to gap filling and trait prediction for macroecology and functional biogeography. Glob. Ecol. Biogeogr. 24, 1510–1521 (2015).Article 

Google Scholar 
Kuppler, J. et al. Global gradients in intraspecific variation in vegetative and floral traits are partially associated with climate and species richness. Glob. Ecol. Biogeogr. 29, 992–1007 (2020).Article 

Google Scholar 
Scheiter, S., Langan, L. & Higgins, S. I. Next-generation dynamic global vegetation models: learning from community ecology. New Phytol. 198, 957–969 (2013).Article 
PubMed 

Google Scholar 
Taubert, F. et al. Confronting an individual-based simulation model with empirical community patterns of grasslands. PLoS ONE 15, e0236546 (2020).Roger, E. & Klistorner, S. (2016) Bioblitzes help science communicators engage local communities in environmental research. J. Sci. Commun. https://doi.org/10.22323/2.15030206 (2016).Legendre, P. & Legendre, L. Numerical Ecology 3rd edn (Elsevier, 2012).Warton, D. I. et al. Smatr 3—an R package for estimation and inference about allometric lines. Methods Ecol Evol 3, 257–259 (2012).Article 

Google Scholar 
Wolf, S. et al. iNaturalist_traits: iNaturalist trait maps version 1 (January 5, 2022) Zenodo https://doi.org/10.5281/zenodo.6671891 (2022). More

De Meyer, M., Robertson, M., Peterson, A. & Mansell, M. Ecological niches and potential geographical distributions of Mediterranean fruit fly (Ceratitis capitata) and Natal fruit fly (Ceratitis rosa). J. Biogeogr. 35, 270–281 (2008).
Google Scholar 
Nguyen, A. D. et al. Trade-offs in cold resistance at the northern range edge of the common woodland ant Aphaenogaster picea (Formicidae). Am. Nat. 194, E151–E163 (2019).Article 
PubMed 

Google Scholar 
Gilioli, G. et al. Non-linear physiological responses to climate change: the case of Ceratitis capitata distribution and abundance in Europe. Biol. Invasions 24, 261–279 (2022).Article 

Google Scholar 
Lancaster, L. T., Dudaniec, R. Y., Hansson, B. & Svensson, E. I. Latitudinal shift in thermal niche breadth results from thermal release during a climate-mediated range expansion. J. Biogeogr. 42, 1953–1963 (2015).Article 

Google Scholar 
Hallas, R., Schiffer, M. & Hoffmann, A. A. Clinal variation in Drosophila serrata for stress resistance and body size. Genet. Res. 79, 141–148 (2002).Article 
PubMed 

Google Scholar 
Hoffmann, A. A., Anderson, A. & Hallas, R. Opposing clines for high and low temperature resistance in Drosophila melanogaster. Ecol. Lett. 5, 614–618 (2002).Article 

Google Scholar 
Ragland, G. & Kingsolver, J. Influence of seasonal timing on thermal ecology and thermal reaction norm evolution in Wyeomyia smithii. J. Evol. Biol. 20, 2144–2153 (2007).Article 
CAS 
PubMed 

Google Scholar 
MacMillan, H. A. & Sinclair, B. J. Mechanisms underlying insect chill-coma. J. Insect Physiol. 57, 12–20 (2011).Article 
CAS 
PubMed 

Google Scholar 
Neilson, E. W. et al. There’sa storm a-coming: Ecological resilience and resistance to extreme weather events. Ecol. Evol. 10, 12147–12156 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
Overgaard, J., Hoffmann, A. A. & Kristensen, T. N. Assessing population and environmental effects on thermal resistance in Drosophila melanogaster using ecologically relevant assays. J. Therm. Biol. 36, 409–416 (2011).Article 

Google Scholar 
Maysov, A. Chill coma temperatures appear similar along a latitudinal gradient, in contrast to divergent chill coma recovery times, in two widespread ant species. J. Exp. Biol. 217, 2650–2658 (2014).Article 
PubMed 

Google Scholar 
David, R. J. et al. Cold stress tolerance in Drosophila: analysis of chill coma recovery in D. melanogaster. J. therm. biol. 23, 291–299 (1998).Article 

Google Scholar 
Overgaard, J. & MacMillan, H. A. The integrative physiology of insect chill tolerance. Annu. Rev. Physiol. 79, 187–208 (2017).Article 
CAS 
PubMed 

Google Scholar 
Andersen, M. K. & Overgaard, J. The central nervous system and muscular system play different roles for chill coma onset and recovery in insects. Comp. Biochem. Physiol. A: Mol. Integr. Physiol. 233, 10–16 (2019).Article 
CAS 

Google Scholar 
Macdonald, S., Rako, L., Batterham, P. & Hoffmann, A. Dissecting chill coma recovery as a measure of cold resistance: evidence for a biphasic response in Drosophila melanogaster. J. Insect Physiol. 50, 695–700 (2004).Article 
CAS 
PubMed 

Google Scholar 
Gibert, P., Moreteau, B., Pétavy, G., Karan, D. & David, J. R. Chill-coma tolerance, a major climatic adaptation among Drosophila species. Evolution 55, 1063–1068 (2001).Article 
CAS 
PubMed 

Google Scholar 
Ayrinhac, A. et al. Cold adaptation in geographical populations of Drosophila melanogaster: phenotypic plasticity is more important than genetic variability. Funct. Ecol. 18, 700–706 (2004).Article 

Google Scholar 
Castañeda, L. E., Lardies, M. A. & Bozinovic, F. Interpopulational variation in recovery time from chill coma along a geographic gradient: a study in the common woodlouse, Porcellio laevis. J. Insect Physiol. 51, 1346–1351 (2005).Article 
PubMed 

Google Scholar 
Tonione, M. A., Cho, S. M., Richmond, G., Irian, C. & Tsutsui, N. D. Intraspecific variation in thermal acclimation and tolerance between populations of the winter ant Prenolepis imparis. Ecol. Evol. 10, 4749–4761 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
Karl, I., Janowitz, S. A. & Fischer, K. Altitudinal life-history variation and thermal adaptation in the copper butterfly Lycaena tityrus. Oikos 117, 778–788 (2008).Article 

Google Scholar 
Ghalambor, C. K., Huey, R. B., Martin, P. R., Tewksbury, J. J. & Wang, G. Are mountain passes higher in the tropics? Janzen’s hypothesis revisited. Integr. Comp. Biol. 46, 5–17 (2006).Article 
PubMed 

Google Scholar 
Addo-Bediako, A., Chown, S. L. & Gaston, K. J. Thermal tolerance, climatic variability and latitude. In Proceedings of the Royal Society of London. Series B: Biological Sciences 267, 739–745 (2000).Poikela, N., Tyukmaeva, V., Hoikkala, A. & Kankare, M. Multiple paths to cold tolerance: the role of environmental cues, morphological traits and the circadian clock gene vrille. BMC ecol. Evol. 21, 1–20 (2021).
Google Scholar 
Andersen, J. L. et al. How to assess Drosophila cold tolerance: chill coma temperature and lower lethal temperature are the best predictors of cold distribution limits. Funct. Ecol. 29, 55–65 (2015).Article 

Google Scholar 
Papadopoulos, N., Katsoyannos, B., Carey, J. & Kouloussis, N. Seasonal and annual occurrence of the Mediterranean fruit fly (Diptera: Tephritidae) in northern Greece. Ann. Entomol. Soc. Am. 94, 41–50 (2001).Article 

Google Scholar 
Malacrida, A. et al. Globalization and fruitfly invasion and expansion: the medfly paradigm. Genetica 131, 1–9 (2007).Article 
CAS 
PubMed 

Google Scholar 
Egartner, A., Lethmayer, C., Gottsberger, R. A. & Blümel, S. In Joint Meeting of the IOBC-WPRS Working Groups “Pheromones and other semiochemicals in integrated production” & “Integrated Protection of Fruit Crops” at. 143–152.Nyamukondiwa, C., Kleynhans, E. & Terblanche, J. S. Phenotypic plasticity of thermal tolerance contributes to the invasion potential of mediterranean fruit flies (Ceratitis capitata). Ecol. Entomol. 35, 565–575 (2010).Article 

Google Scholar 
Weldon, C. W., Terblanche, J. S. & Chown, S. L. Time-course for attainment and reversal of acclimation to constant temperature in two Ceratitis species. J. Therm. Biol. 36, 479–485 (2011).Article 

Google Scholar 
Pujol-Lereis, L. M., Rabossi, A. & Quesada-Allué, L. A. Analysis of survival, gene expression and behavior following chill-coma in the medfly Ceratitis capitata: effects of population heterogeneity and age. J. Insect Physiol. 71, 156–163 (2014).Article 
CAS 
PubMed 

Google Scholar 
Pujol-Lereis, L. M., Fagali, N. S., Rabossi, A., Catalá, Á. & Quesada-Allué, L. A. Chill-coma recovery time, age and sex determine lipid profiles in Ceratitis capitata tissues. J. Insect Physiol. 87, 53–62 (2016).Article 
CAS 
PubMed 

Google Scholar 
Weldon, C. W., Nyamukondiwa, C., Karsten, M., Chown, S. L. & Terblanche, J. S. Geographic variation and plasticity in climate stress resistance among southern African populations of Ceratitis capitata (Wiedemann)(Diptera: Tephritidae). Sci. Rep. 8, 1–13 (2018).Article 
CAS 

Google Scholar 
Nyamukondiwa, C., Weldon, C. W., Chown, S. L., le Roux, P. C. & Terblanche, J. S. Thermal biology, population fluctuations and implications of temperature extremes for the management of two globally significant insect pests. J. Insect Physiol. 59, 1199–1211 (2013).Article 
CAS 
PubMed 

Google Scholar 
Mitchell, K. A., Boardman, L., Clusella-Trullas, S. & Terblanche, J. S. Effects of nutrient and water restriction on thermal tolerance: A test of mechanisms and hypotheses. Comp. Biochem. Physiol. A: Mol. Integr. Physiol. 212, 15–23 (2017).Article 
CAS 

Google Scholar 
Hoffmann, A. A. & Ross, P. A. Rates and patterns of laboratory adaptation in (mostly) insects. J. Econ. Entomol. 111, 501–509 (2018).Article 
PubMed 

Google Scholar 
Popa-Báez, Á. -D. et al. Climate stress resistance in male Queensland fruit fly varies among populations of diverse geographic origins and changes during domestication. BMC Genet. 21, 1–19 (2020).Article 

Google Scholar 
Beck, H. E. et al. Present and future Köppen-Geiger climate classification maps at 1-km resolution. Sci. Data 5, 1–12 (2018).Article 

Google Scholar 
Kozak, K. H., Graham, C. H. & Wiens, J. J. Integrating GIS-based environmental data into evolutionary biology. Trends Ecol. Evol. 23, 141–148 (2008).Article 
PubMed 

Google Scholar 
Oyen, K. J. et al. Body mass and sex, not local climate, drive differences in chill coma recovery times in common garden reared bumble bees. J. Comp. Physiol. B. 191, 843–854 (2021).Article 
PubMed 

Google Scholar 
Angert, A. L., Bontrager, M. G. & Ågren, J. What do we really know about adaptation at range edges?. Annu. Rev. Ecol. Evol. Syst. 51, 341–361 (2020).Article 

Google Scholar 
Terblanche, J. S. & Hoffmann, A. A. Validating measurements of acclimation for climate change adaptation. Curr. Opin. insect sci. 41, 7–16 (2020).Article 
PubMed 

Google Scholar 
Kourti, A. Patterns of variation within and between Greek populations of Ceratitis capitata suggest extensive gene flow and latitudinal clines. J. Econ. Entomol. 97, 1186–1190 (2004).Article 
CAS 
PubMed 

Google Scholar 
Hangartner, S., Lasne, C., Sgrò, C. M., Connallon, T. & Monro, K. Genetic covariances promote climatic adaptation in Australian Drosophila. Evolution 74, 326–337 (2020).Article 
PubMed 

Google Scholar 
Bontrager, M. & Angert, A. L. Gene flow improves fitness at a range edge under climate change. Evol. Let. 3, 55–68 (2019).Article 

Google Scholar 
Liu, Q. et al. Extension of the growing season increases vegetation exposure to frost. Nat. Commun. 9, 1–8 (2018).
Google Scholar 
Schwartz, M. D., Ahas, R. & Aasa, A. Onset of spring starting earlier across the Northern Hemisphere. Glob. Change Biol. 12, 343–351 (2006).Article 

Google Scholar 
Ma, Q., Huang, J. G., Hänninen, H. & Berninger, F. Divergent trends in the risk of spring frost damage to trees in Europe with recent warming. Glob. Change Biol. 25, 351–360 (2019).Article 

Google Scholar 
Unterberger, C. et al. Spring frost risk for regional apple production under a warmer climate. PLoS ONE 13, e0200201 (2018).Article 
MathSciNet 
PubMed 
PubMed Central 

Google Scholar 
Manrakhan, A., Daneel, J.-H., Stephen, P. R. & Hattingh, V. Cold Tolerance of Immature Stages of Ceratitis capitata and Bactrocera dorsalis (Diptera: Tephritidae). J. Econ. Entomol. 115(2), 482–492 (2022).Article 
PubMed 

Google Scholar 
Papadopoulos, N. T., Carey, J. R., Katsoyannos, B. I. & Kouloussis, N. A. Overwintering of the mediterranean fruit fly (Diptera: Tephritidae) in Northern Greece. Ann. Entomol. Soc. Am. 89, 526–534 (1996).Article 

Google Scholar 
Papadopoulos, N. T., Katsoyannos, B. I. & Carey, J. R. Temporal changes in the composition of the overwintering larval population of the Mediterranean fruit fly (Diptera: Tephritidae) in Northern Greece. Ann. Entomol. Soc. Am. 91, 430–434 (1998).Article 

Google Scholar 
Katsoyannos, B. I., Kouloussis, N. A. & Carey, J. R. Seasonal and annual occurrence of Mediterranean fruit flies (Diptera: Tephritidae) on Chios Island, Greece: Differences between two neighboring citrus orchards. Ann. Entomol. Soc. Am. 91, 43–51 (1998).Article 

Google Scholar 
Mavrikakis, P. G., Economopoulos, A. P. & Carey, J. R. Continuous winter reproduction and growth of the mediterranean fruit fly (Diptera: Tephritidae) in Heraklion, crete Southern Greece. Environ. Entomol. 29, 1180–1187 (2000).Article 

Google Scholar 
Israely, N., Ziv, Y. & Oman, S. D. Spatiotemporal distribution patterns of Mediterranean fruit fly (Diptera: Tephritidae) in the central region of Israel. Ann. Entomol. Soc. Am. 98, 77–84 (2005).Article 

Google Scholar 
Bahrndorff, S., Lauritzen, J. M., Sørensen, M. H., Noer, N. K. & Kristensen, T. N. Responses of terrestrial polar arthropods to high and increasing temperatures. J. Exp. Biol. 224, jeb230797 (2021).Article 
PubMed 

Google Scholar 
Sinclair, B. J. & Roberts, S. P. Acclimation, shock and hardening in the cold. J. Therm. Biol. 30, 557–562 (2005).Article 

Google Scholar 
Bahrndorff, S., Gertsen, S., Pertoldi, C. & Kristensen, T. N. Investigating thermal acclimation effects before and after a cold shock in Drosophila melanogaster using behavioural assays. Biol. J. Lin. Soc. 117, 241–251 (2016).Article 

Google Scholar 
Sarmad, M., Ishfaq, A., Arif, H. & Zaka, S. M. Effect of short-term cold temperature stress on development, survival and reproduction of Dysdercus koenigii (Hemiptera: Pyrrhocoridae). Cryobiology 92, 47–52 (2020).Article 
PubMed 

Google Scholar 
Steyn, V. M., Mitchell, K. A., Nyamukondiwa, C. & Terblanche, J. S. Understanding costs and benefits of thermal plasticity for pest management: Insights from the integration of laboratory, semi-field and field assessments of Ceratitis capitata (Diptera: Tephritidae). Bull. Entomol. Res., 1–11 (2022).Davis, H. E., Cheslock, A. & MacMillan, H. A. Chill coma onset and recovery fail to reveal true variation in thermal performance among populations of Drosophila melanogaster. Sci. Rep. 11, 1–10 (2021).Article 

Google Scholar 
Noh, S., Everman, E. R., Berger, C. M. & Morgan, T. J. Seasonal variation in basal and plastic cold tolerance: Adaptation is influenced by both long-and short-term phenotypic plasticity. Ecol. Evol. 7, 5248–5257 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
Bruins, H. J. Ancient desert agriculture in the Negev and climate-zone boundary changes during average, wet and drought years. J. Arid Environ. 86, 28–42 (2012).Article 

Google Scholar 
Hoffmann, A. A., Sørensen, J. G. & Loeschcke, V. Adaptation of Drosophila to temperature extremes: Bringing together quantitative and molecular approaches. J. Therm. Biol. 28, 175–216 (2003).Article 

Google Scholar 
Kawecki, T. J. & Ebert, D. Conceptual issues in local adaptation. Ecol. Lett. 7, 1225–1241 (2004).Article 

Google Scholar 
Nyamukondiwa, C. & Terblanche, J. S. Thermal tolerance in adult mediterranean and Natal fruit flies (Ceratitis capitata and Ceratitis rosa): Effects of age, gender and feeding status. J. Therm. Biol. 34, 406–414 (2009).Article 

Google Scholar 
Team, R. C. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria (2021).Fick, S. E. & Hijmans, R. J. WorldClim 2: New 1-km spatial resolution climate surfaces for global land areas. Int. J. Climatol. 37, 4302–4315 (2017).Article 

Google Scholar 
Mazerolle, M. J. Model selection and multimodel inference using the AICcmodavg package (2020).Therneau, T. A Package for Survival Analysis in R. R Package Version 3.2-13.(2021. (2021).Kassambara, A., Kosinski, M., Biecek, P. & Fabian, S. Survminer: Drawing Survival Curves using’ggplot2′. R package version 0.4. 9. 2021. (2021).Lenth, R. V. Emmeans: Estimated Marginal Means, aka Least-Squares Means. R package version 1.7.2. (2022). More

Study areaThe present study was conducted in Uttarkashi district, Uttarakhand, located between 38° 28′ to 31°28′ N latitude and 77°49′ to 79°25′ E longitude with an area of about 8016 km2, covering primarily hilly terrain with an altitudinal range of 715–6717 m (Fig. 3). The terrain is mountainous, consisting of undulating hill ranges and narrow valleys with temperate climatic conditions. The district lies in the upper catchment of two major rivers of India, viz., the Ganges (Bhagirathi towards upstream) and the Yamuna. The major vegetation types of the study area are Himalayan moist temperate forest, sub-alpine forest and alpine scrub59. The Uttarkashi district forests are managed under three Forest Divisions viz., (i) Uttarkashi Forest Division (ii) Upper Yamuna Badkot Forest Division and (iii) Tons Forest Division) with two Protected Areas (PAs) (i) Gangotri National Park and (ii) Govind Pashu Vihar National Park. The forested habitats of the study landscape are home to top conservation priority species, including Asiatic Black bear (Ursus thibetanus), Musk deer (Moschus spp.), Common leopard (Panthera pardus), Himalayan brown bear (Ursus arctos isabellinus) and Western Tragopan (Tragopan melanocephalus), Himalayan monal (Lophophorus impejanus). The study was conducted after a study permit issued by the Chief Wildlife Warden, Forest Department, Uttarakhand government, vide letter no. 848/5-6 dated 31/08/2019, we have not handled the species for doing research. Instead, remote camera traps have been used for collecting the data with the permission of the Chief Wildlife Warden, Government of Uttarakhand. Further, informed consent was taken before interviewing the local communities. The data was collected according to the institutional guidelines and approved by the Research Advisory and Monitoring Committee of the Zoological Survey of India.Figure 3Map of the study area Uttarkashi, Uttarakhand. ArcGIS 10.6 (ESRI, Redlands, CA) was used to create the map. (Map created using ArcGIS 10.6; http://www.esri.com).Full size imageSampling protocolThe basic sampling protocol and assumptions for multi-species occupancy modelling are identical to the single-species case7. Briefly, a set of 62 intensive sites, were randomly selected, and each site i was surveyed j times. During each survey, detection/non-detection of S focal species was recorded. Additionally, direct or indirect evidences of species presence from the different areas were also recorded.Data collectionThe complete study area was divided into 10 × 10 km grids, consisting of n = 60 grids. Based on the reconnaissance survey, out of these 60 grids, we selected 25 girds that were accessible to conduct the survey and have the species presence. Further, these grids were divided into 2 × 2 km grids to maximize our effort so that all logistically accessible grids could be covered, and we conducted intensive sampling in N = 62 grids after excluding the grids with human settlements. T The field surveys were conducted during 2018–2019, and a team of researchers systematically visited selected grids to collect data on the detection/non-detection of these ungulates. A total of 62 camera traps were deployed in selected grids, and 650 km were traversed, accounting for N = 54 trails in these sampled grids. These camera traps were visited once in every fifteen days for replacing the batteries as well as documenting the presence of the species through the sign surveys. The ultra-compact SPYPOINT FORCE-11D trail camera (SPYPOINT, GG Telecom, Canada, QC) and Browning trail camera (Defender 850, 20 MP, Prometheus Group, LLC Birmingham, Alabama, https://browningtrailcameras.com) camera traps were used to detect the presence/absence of ungulate species. The cameras were mounted 40–60 cm above ground on natural trails without lures.Data explorationWhile deploying camera traps, we also noted habitat variables through on-site observation such as distance to the village and human disturbance. We tested site covariates for collinearity and discarded one of a pair if the Pearson’s correlation was greater than 0.760. Hence, we assumed each of the site covariates could influence the occupancy and detectability of these ungulates.CovariatesWe hypothesized that habitat variables may influence these ungulates’ occupancy and detection probability. A total of 21 variables were extracted either from the field or using the ArcGIS v. 10.6 software (ESRI, Redlands, CA), and only 14 were retained after collinearity testing60 (Table 3). These covariates were classified into the following categories (Topographic variables, Habitat variables and anthropogenic variables). The topographic variables (elevation, slope and aspect) were generated using 30× resolution SRTM (Shuttle Radar Topography Mission) image downloaded from EarthExplorer (https://earthexplorer.usgs.gov/). The habitat/ land cover classification was carried out using Landsat 8 satellite imagery (Spatial resolution = 30 m) downloaded from Global Land Cover Facility by following the methodology suggested by61 using the ArcGIS v. 10.6 software (ESRI, Redlands, CA). The study area was classified into nine Land use/land cover (LULC) classes viz., West Himalayan Sub-alpine birch/fir Forest (FT 188), West Himalayan upper oak/fir forest (FT 162), West Himalayan Dry juniper forest (FT 180), Ban oak forest (FT 152), Moist Deodar Forest (FT 155), Western mixed coniferous forest (FT 156), Moist temperate Deciduous Forest (FT 157) which were used for further analysis considering their importance to species ecology and behavior60. The values for all the covariates were extracted at 30 m resolution, and a single value per site was obtained by averaging all the pixel values within each sampling site (camera trap locations).Table 3 Habitat variables used for multi species occupancy analysis of three ungulate species in Uttarkashi, Uttarakhand.Full size tableOccupancy modelling frameworkWe used multi-species occupancy modelling62 of barking deer, goral and sambar to estimate the probability of the species (s) occurred within the area (i) sampled during our survey period (j), for accounting the imperfect detection of the species8. Distinguishing the true presence/absence of a species from detection/non-detection (i.e., species present and captured or species present but not captured) requires spatially or temporally replicated data. We used camera stations to record the presence/absence of species along with sign survey in all the studied grids. The camera traps were placed along the trail/transects in the studied grids hence each grid needs to be visited once in every fifteen days to check the camera traps as well as to document the presence of the studied species. Therefore, we treated 15 trap nights as one sampling occasion at a particular camera station resulting in ~ 7 sampling occasions per camera station.Our aim was to record the presence/ absence of the species at a particular gird hence we incorporated sign survey data if the species was not detected in camera station but recorded through sign survey. We pooled the presence/absence data in a single sheet of each species following6 and fitted occupancy and detectability models using programme Mark63,64. We model the species (s) presence (ysij = 1) and absence (ysij = 0) at site i during survey j, and the sampling protocol was identical to single species case65, where the Bernoulli random variable was conditional on the presence of species s (Zs = 1) following6$${text{y}}_{sij} sim {text{ Bernoulli}}left( {{text{p}}_{sij} {text{z}}_{si} } right),$$
where Psij represents the probability of detecting species S during replicate survey j at site i and Zsi = presence or absence of species s at site i.Furthermore, we model the latent occupancy state of species s at site i as a multivariate Bernoulli random variable:$${text{Z}}_{i} sim {text{MVB}}left( {uppsi _{i} } right)$$
where Zi = {Z1i, Z2i….., ZSi} is an S-dimensional vector of 1’s and 0’s denoting the latent occupancy state of all S species and (ψi) is a 2S-dimensional vector denoting the probability of all possible sequences of 1’s and 0’s Zi can attain such that ∑ ψi = 1 with corresponding probability mass function (PMF) adopted from6,64.$$fleft( {{text{Z}}_{i} } right) = {text{ exp}}left( {left( {{text{Z}}_{i} {text{log}}(uppsi_{{text{i}}} {1}/uppsi_{{text{i}}} 0} right) , + {text{ log}}left( {uppsi_{{text{i}}} 0} right)} right).$$The quantity f = log (ψi1/ψi0), is the log odds species S occupies a site often referred to as a ‘natural parameter’.Since we are modeling three ungulate species (S = 3), 2S = 23 the possible encounter histories included in the dataset were eight, if neither of the two species were detected the value of ‘00’ was assigned; similarly ‘01’ indicates detection of species 1; ‘02’ indicates detection of species 2; ‘03’ indicates detection of both the species; ‘04’ indicates detection of species 3; ‘05’ indicates detection of species 1 and species 3; ‘06’ indicates detection of species 2 and species 3 and ‘07’ indicates detection of all the three species. We modelled constant occupancy and detection probability for each of the three species. Hence, we specified 6 f and p parameters, an intercept (β) for each of one-way f parameter and detection parameter p following64.$$f_{{1}}=upbeta_{{{1},}} ;;{text{p}}=upbeta_{{4}}$$$$f_{{2}} = upbeta_{{{2},}} ;{text{p }} = , upbeta 5$$$$f_{{3}} = , upbeta_{{{3},}}; {text{p }} = , upbeta_{{6}}$$We fit a set of models including the detection probability as a constant, p(.), and variable function to occupancy ψ(covariate) for site-specific covariates and models include occupancy as constant ψ(.) and variable function of the detection p(covariates) for the respective site covariates.As we have assumed the independence among all three species, the model shows marginal occupancy probabilities of species 1, species 2 and species 3 varies as a function of environmental variables. We incorporated site-level characteristics affecting species-specific occurrence (f1: occupancy of species 1, f2: occupancy of species 2, & f3: occupancy of species 3) and detection probabilities using a generalized linear modelling approach42. This requires 9 parameters: an intercept (β1, β3, β5) and slope (β2, β4, β6) coefficient for each 1-way f parameter f1, f2, f3 and an intercept parameter for each detection parameter (β7, β8, β9). Below mentioned is the model for 1-way f parameters.$$f_{{1}} = , upbeta_{{{1 } + }} upbeta_{{2}} left( {{text{Covariate}}} right),;;{text{ p }} = , upbeta_{{7}}$$$$f_{{2}} = , upbeta_{{{3 } + }} upbeta_{{4}} left( {{text{Covariate}}} right),;;{text{ p}} = , upbeta_{{8}}$$$$f_{{3}} = , upbeta_{{5}} + , upbeta_{{6}} left( {{text{Covariate}}} right),;;{text{ p }} = , upbeta_{{9}} .$$All covariates were standardized before model fitting. We fitted the most complex model to each species and considered all possible combinations of covariates using the logit link function. Our rationale for including these variables in the occupancy and detectability component of the model was that we expected these variables to influence the occupancy and detectability of the study species.Since multi-species occupancy simultaneously model environmental variables, & interspecific interactions. Further it also allows to understand the influence of environmental variables on one species occupancy, in the presence or absence of other sympatric species64. Hence, we also modeled two species occur together as a function of covariates. We examined how the variables of each camera site influenced the pair-wise interaction of the three ungulate species. This model assumes that the conditional probability of one species varies in the presence or absence of other species. We assumed f123: co-occurrence of species 1, species 2 & species 3 = 0, hence we did not include higher-order interactions in any of our models, we assumed the conditional probability of 3 species occurred together was purely a function of species-specific (f1, f2, f3) and pair-wise interaction (f12: co-occurrence of species1 & species 2, f13: co-occurrence of species 1 & species 3, f23: co-occurrence of species 2 & species 3) parameters. We modeled pair-wise interaction of species varies as a function of environmental variables keeping detection probability constant. Hence, we specified 15 f and p parameters, an intercept and slope coefficient for each of the one-way (f1, f2, f3) and the two-way f parameters (f12, f13, and f23); as well as an intercept parameter for each of the detection models. The model equation below implies for 2-way f parameters:$$f_{{{12}}} = , upbeta_{{{7 } + }} upbeta_{{8}} left( {{text{Covariate}}} right),;;{text{ p }} = , upbeta_{{{13}}}$$$$f_{{{13}}} = , upbeta_{{{9 } + }} upbeta_{{{1}0}} left( {{text{Covariate}}} right),;;{text{ p }} = , upbeta_{{{14}}}$$$$f_{{{23}}} = , upbeta_{{{11 } + }} upbeta_{{{12}}} left( {{text{Covariate}}} right),;;{text{ p }} = , upbeta_{{{15}}} .$$We also fitted models including co-occurrence and detection probability of a species varies as a function of environmental variables. Hence, we specified 18 f and p parameters, an intercept and slope coefficient for each of one-way (f1, f2, f3) and two-way f parameters (f12, f13, f23); and an intercept as well as the slope parameters for each of the detection models. The model equation below implies for 2-way f parameters:$$f_{{{12}}} = , upbeta_{{{7 } + }} upbeta_{{8}} left( {{text{Covariate}}} right),{text{ p }} = , upbeta_{{{13 } + }} upbeta_{{{14}}} left( {{text{covariate}}} right)$$$$f_{{{13}}} = , upbeta_{{{9 } + }} upbeta_{{{1}0}} left( {{text{Covariate}}} right),{text{ p }} = , upbeta_{{{15}}} + , upbeta_{{{16}}} left( {{text{covariate}}} right)$$$$f_{{{23}}} = , upbeta_{{{11 } + }} upbeta_{{{12}}} left( {{text{Covariate}}} right),{text{ p }} = , upbeta_{{{17}}} + , upbeta_{{{18}}} left( {{text{covariate}}} right)$$A total of 38 models were run to test the influence of environmental variables on occupancy and detection probability of species-specific (f1, f2, f3) and pair-wise interaction of the three ungulate species. The best-supported model was identified by selecting the model with the lowest AICc value and highest model weights66, where higher model weights indicate a better fit of the model to the data. Second-Order Information Criterion (AICc)67 values were used to rank the occupancy models, and all the models whose ΔAICc  More

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.This is a summary of: Marin, F. R. et al. Protecting the Amazon forest and reducing global warming via agricultural intensification. Nat. Sustain. https://doi.org/10.1038/s41893-022-00968-8 (2022). More