Philippa Kaur

Department of Ecology and Evolutionary Biology and Eversource Energy Center, University of Connecticut, Storrs, CT, USACory MerowDepartment of Ecology and Evolutionary Biology, University of Arizona, Tucson, AZ, USABrad Boyle & Brian J. EnquistDepartment of Geography, Florida State University, Tallahassee, FL, USAXiao FengBiodiversity and Biocomplexity Unit, Okinawa Institute of Science and Technology Graduate University, Onna, JapanJamie M. KassDepartment of Geography, University at Buffalo, Buffalo, NY, USABrian S. Maitner & Adam M. WilsonSchool of Biology and Ecology, University of Maine, Orono, ME, USABrian McGillMitchell Center for Sustainability Solutions, University of Maine, Orono, ME, USABrian McGillCenter for Macroecology, Evolution and Climate, Globe Institute, University of Copenhagen, Copenhagen, DenmarkHannah OwensFlorida Museum of Natural History, University of Florida, Gainesville, FL, USAHannah OwensDepartment of Biological Sciences, Purdue University, West Lafayette, IN, USADaniel S. ParkPurdue Center for Plant Biology, Purdue University, West Lafayette, IN, USADaniel S. ParkDepartment of Environmental Systems Science, Institute of Integrative Biology, ETH Zürich, Zurich, SwitzerlandAndrea PazDepartment of Biology, City College of the City University of New York, New York, NY, USAGonzalo E. Pinilla-BuitragoPhD Program in Biology, Graduate Center of the City University of New York, New York, NY, USAGonzalo E. Pinilla-BuitragoDepartment of Ecology and Evolutionary Biology, University of Connecticut, Storrs, CT, USAMark C. UrbanCenter of Biological Risk, University of Connecticut, Storrs, CT, USAMark C. UrbanDepartamento de Ecoloxía e Bioloxía Animal, Universidade de Vigo, Vigo, SpainSara Varela More

Breed differences in animal feed conversion and economic trait performanceThroughout the feed intake measurement period, summary statistics shows animals on test had an average DMI of 1.11 kg/d (SD = 0.18), ADG of 0.27 kg/d (SD = 0.1), FCR of 4.04 kg of DMI/ Kg of ADG (SD = 0.1), start weight of 29.60 kg (SD = 3.7), final live weight of 46.00 kg (SD = 2.9), carcass weight of 20.20 kg (SD = 1.6), and a KO% of 44.1% (SD = 2.3). Average daily gain (P = 0.005), FCR (P = 0.035), CW (P  More

Finkelman, S., Navarro, S., Rindner, M. & Dias, R. Effect of low pressure on the survival of Trogoderma granarium Everts, Lasioderma serricorne (F.) and Oryzaephilus surinamensis (L.) at 30°C. J. Stored. Prod. Res. 42, 23–30 (2006).Article 

Google Scholar 
Hosseininaveh, V., Bandani, A., Azmayeshfard, P., Hosseinkhani, S. & Kazzazi, M. Digestive proteolytic and amylolytic activities in Trogoderma granarium Everts (Dermestidae: Coleoptera). J. Stored. Prod. Res. 43, 515–522 (2007).Article 
CAS 

Google Scholar 
Burges, H. D. Development of the khapra beetle, Trogoderma granarium, in the lower part of its temperature range. J. Stored. Prod. Res. 44, 32–35 (2008).Article 

Google Scholar 
Hagstrum, D. W. & Subramanyam, B. Stored-Product Insect Resource 1–518 (AACC International Inc, 2009).Book 

Google Scholar 
Beal, R. S. Synopsis of the economic species of Trogoderma occurring in the United States with description of a new species (Coleoptera: Dermestidae). Ann. Entomol. Soc. Am. 49, 559–566 (1956).Article 

Google Scholar 
Day, C. & White, B. Khapra beetle, Trogoderma granarium interceptions and eradications in Australia and around the world. Crawley, School of Agricultural and Resource Economics, University of Western Australia, SARE Working paper 1609, (2016).Kerr, J. A. Khapra beetle returns. Pest Control 49, 24–25 (1981).
Google Scholar 
Stibick, J.N. New pest response guidelines: khapra beetle. US Department of Agriculture, Marketing and Regulatory Programs, Animal and Plant Health Inspection Service, Riverdale, pp. 114 (2009).Myers, S. W. & Hagstrum, D. W. Quarantine. In Stored Product Protection (eds Hagstrum, D. W. et al.) 297–304 (Kansas State University Agricultural Experiment Station and Cooperative Extension Service, 2012).
Google Scholar 
Athanassiou, C. G., Phillips, T. W. & Wakil, W. Biology and control of the khapra beetle, Trogoderma granarium, a major quarantine threat to global food security. Annu. Rev. Entomol. 64, 131–148 (2019).Article 
CAS 
PubMed 

Google Scholar 
Barak, A. V. Development of a new trap to detect and monitor khapra beetle (Coleoptera: Dermestidae). J. Econ. Entom. 82, 1470–1477 (1989).Article 

Google Scholar 
Gerken, A. R. & Campbell, J. F. Life history changes in Trogoderma variabile and T. inclusum due to mating delay with implications for mating disruption as a management tactic. Ecol. Evol. 8, 2428–2439 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Partida, G. J. & Strong, R. G. Comparative studies on the biologies of six species of Trogoderma: T variabile. Ann. Entomol. Soc. Am. 68, 115–125 (1975).Article 

Google Scholar 
Strong, R. G. Comparative studies on the biologies of six species of Trogoderma: T inclusum. Ann. Entomol. Soc. Am. 68, 91–104 (1975).Article 

Google Scholar 
Phillips, T. W., Pfannenstiel, L. & Hagstrum, D. Survey of Trogoderma species (Coleoptera: Dermestidae) associated with international trade of dried distiller’s grains and solubles in the USA. Julius Kühn Archiv. 463, 233–238 (2008).
Google Scholar 
Hadaway, A. The biology of the beetles, Trogoderma granarium Everts and Trogoderma versicolor (Creutz). Bull. Entomol. Res. 46, 781–796 (1956).Article 
CAS 

Google Scholar 
Phillips, T.W., Pfannenstiel, L. & Hagstrum, D. Survey of Trogoderma species (Coleoptera: Dermestidae) associated with international trade of dried distiller’s grains and solubles in the USA. In: Adler CS, Opit G, Fürstenau B, Müller-Blenkle C, Kern P, Arthur FH et al., editors. Proceedings of the 12th International Working Conference on Stored Product Protection; Vol. 1, Quedlinburg, Julius-Kühn-Archiv, pp. 233–238 (2018).Gorham, J.R. Insect and Mite Pests in Food: An Illustrated Key. Vol. 1 and 2. US Department of Agriculture, Agricultural Research Service (1991).Olson, R. L. O., Farris, R. E., Barr, N. B. & Cognato, A. I. Molecular identification of Trogoderma granarium (Coleoptera: Dermestidae) using the 16S gene. J. Pest Sci. 87, 701–710 (2014).Article 

Google Scholar 
Furui, S., Miyanoshita, A., Imamura, T., Minegishi, Y. & Kokutani, R. Qualitative real-time PCR identification of the khapra beetle, Trogoderma granarium (Coleoptera: Dermestidae). Appl. Entomol. Zool. 54, 101–107 (2019).Article 
CAS 

Google Scholar 
Rako, L. et al. A LAMP (loop-mediated isothermal amplification) test for rapid identification of Khapra beetle (Trogoderma granarium). Pest Manag. Sci. 77, 5509–5521 (2021).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Castañé, C., Agustí, N., del Estal, P. & Riudavets, J. Survey of Trogoderma spp. in Spanish mills and warehouses. J. Stored. Prod. Res. 88, 101661 (2020).Article 

Google Scholar 
Trujillo-González, et al. Detection of khapra beetle environmental DNA using portable technologies in Australian biosecurity. Front. Insect Sci. 2, e795379 (2022).Article 

Google Scholar 
Svec, D., Tichopad, A., Novosadova, V., Pfaffl, M. W. & Kubista, M. How good is a PCR efficiency estimate: Recommendations for precise and robust qPCR efficiency assessments. Biomol. Detect. Quantif. 3, 9–16 (2015).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Taylor, S. C. et al. The Ultimate qPCR experiment: Producing publication quality, reproducible data the first time. Trends Biotechnol. 37, 761–774 (2019).Article 
CAS 
PubMed 

Google Scholar 
Van Holm, W. et al. A viability quantitative PCR dilemma: Are longer amplicons better?. Appl. Environ. Microbiol. 87, e0265320 (2021).Article 
PubMed 

Google Scholar 
Ratnasingham, S. & Hebert, P. D. N. BOLD: The barcode of life data system (wwwbarcodinglifeorg). Mol. Ecol. Notes 7, 355–364 (2007).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Wittwer, C. T. & Kusakawa, N. Real-time PCR. In Molecular microbiology: Diagnostic principles and practice (eds Persing, D. H. et al.) 71–84 (ASM Press, 2004).
Google Scholar 
Stewart, D. et al. A needle in a haystack: A multigene TaqMan assay for the detection of Asian gypsy moths in bulk pheromone trap samples. Biol. Invasions 21, 1843–1856 (2019).Article 

Google Scholar 
Butterwort, V. et al. A DNA extraction method for insects from sticky traps: Targeting a low abundance pest, Phthorimaea absoluta (Lepidoptera: Gelechiidae), in mixed species communities. J. Econ. Entom. 115, 844–851 (2022).Article 
CAS 
PubMed 

Google Scholar 
Carew, M. E., Coleman, R. A. & Hoffmann, A. A. Can non-destructive DNA extraction of bulk invertebrate samples be used for metabarcoding?. PeerJ 6, e4980 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Domingue, M.J. et al. Outcome of direct competition between Trogoderma granarium and Trogoderma inclusum over varying commodities, temperatures, and experimental duration. In Submission to Scientific Reports.Zieritz, A. et al. Development and evaluation of hotshot protocols for cost- and time-effective extraction of PCR-ready DNA from single freshwater mussel larvae (Bivalvia: Unionida). J. Molluscan Stud. 84, 198–201 (2018).Article 

Google Scholar 
Djoumad, A. et al. Development of a qPCR-based method for counting overwintering spruce budworm (Choristoneura fumiferana) larvae collected during fall surveys and for assessing their natural enemy load: A proof-of-concept study. Pest Manag. Sci. 78, 336–343 (2022).Article 
CAS 
PubMed 

Google Scholar 
Chen, H., Rangasamy, M., Tan, S. Y., Wang, H. & Siegfried, B. D. Evaluation of five methods for total DNA extraction from western corn rootworm beetles. PLoS ONE 5, e11963 (2010).Article 
ADS 
PubMed 
PubMed Central 

Google Scholar 
Beckmann, J. S. & Soller, M. Restriction fragment length polymorphisms in genetic improvement: Methodologies, mapping and costs. Theor. Appl. Genet. 67, 35–43 (1983).Article 
CAS 
PubMed 

Google Scholar 
Arimoto, M., Satoh, M., Uesugi, R. & Osakabe, M. PCR-RFLP analysis for identification of Tetranychus spider mite species (Acari: Tetranychidae). J. Econ. Entom. 106, 661–668 (2013).Article 
CAS 
PubMed 

Google Scholar 
Vezenegho, S. B. et al. Discrimination of 15 Amazonian anopheline mosquito species by polymerase chain reaction—Restriction fragment length polymorphism. J. Med. Entomol. 59, 1060–1064 (2022).Article 
CAS 
PubMed 

Google Scholar 
Beal, R. S. Annotated checklist of Nearctic Dermestidae with revised key to the genera. Coleopt. Bull. 57, 391–404 (2003).Article 

Google Scholar 
Folmer, O., Black, M., Hoeh, W., Lutz, R. & Vrijenhoek, R. DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates. Mol. Mar. Biol. Biotechnol. 3, 294–299 (1994).CAS 
PubMed 

Google Scholar 
Liu, H. & Mottern, J. An old remedy for a new problem? Identification of Ooencyrtus kuvanae (Hymenoptera: Encyrtidae), an egg parasitoid of Lycorma delicatula (Hemiptera: Fulgoridae) in North America. J. Insect Sci. 17, 18 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
Simon, C. et al. Evolution, weighing, and phylogenetic utility of mitochondrial gene sequences and a compilation of conserved polymerase chain reaction primers. Ann. Entomol. Soc. Am. 87, 651–701 (1994).Article 
CAS 

Google Scholar 
Dowton, M. & Austin, A. D. Evidence for AT-transversion bias in wasp (Hymenoptera: Symphyta) mitochondrial genes and its implications for the origin of parasitism. J. Mol. Evol. 44, 398–405 (1997).Article 
ADS 
CAS 
PubMed 

Google Scholar 
Untergasser, A. et al. Primer3—New capabilities and interfaces. Nucleic Acids Res. 40, e115 (2012).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Ye, J. et al. Primer-BLAST: A tool to design target-specific primers for polymerase chain reaction. BMC Bioinform. 13, 134 (2012).Article 
CAS 

Google Scholar 
Süss, B., Flekna, G., Wagner, M. & Hein, I. Studying the effect of single mismatches in primer and probe binding regions on amplification curves and quantification in real-time PCR. J. Microbiol. Methods 76, 316–319 (2009).Article 
PubMed 

Google Scholar 
Stadhouders, R. et al. The effect of primer-template mismatches on the detection and quantification of nucleic acids using the 5’ nuclease assay. J. Mol. Diagn. 12, 109–117 (2010).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Stewart, D. et al. A multi-species TaqMan PCR assay for the identification of Asian gypsy moths (Lymantria spp.) and other invasive Lymantriines of biosecurity concern to North America. PLoS ONE 11, e0160878 (2016).Article 
PubMed 
PubMed Central 

Google Scholar 
Bustin, S. A. et al. The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments. Clin. Chem. 55, 1–12 (2009).Article 

Google Scholar  More

Bergmann’s ruleVariation in avian body size has arisen through millions of years of evolution43, and our data reflects this by showing that log body mass is strongly predicted by phylogeny (Supplementary Table 1). Yet, avian body size also shows large geographical variation (Fig. 1a), and our analysis provides strong support for Bergmann’s rule across the global community of birds. Phylogenetic linear mixed models indicated that the temperature variables explain from 9.0% to 11.8% of the variance in log-transformed body size (estimated with r-squared; Fig. 1b). These models are substantially better supported than the null model and the model with latitude alone (Fig. 1b), suggesting that the observed geographical pattern is linked to thermoregulation. All of these temperature models indicate that temperature negatively correlates with body size (Fig. 1c and Supplementary Fig. 1), as predicted by Bergmann’s rule.Fig. 1: Global test of Bergmann’s rule across 9962 (99.7%) avian species.Distribution of log-transformed body mass across species geographic ranges, shown as their geometric centroids (a). Model selection procedure for predicting log body mass (b), with six temperature measures assessed within species geographic ranges, as sole fixed effects; AIC—Akaike Information Criterion, r2—coefficient of determination. An exemplar Bergmann’s model (c), showing decreasing body size with max temperature of all months; see Supplementary Fig. 1 for surrogate models based on the other temperature measures (evaluated in b). The shaded area around the trend line is simple shading to facilitate reading. The p values refer to the significance of temperature effect and whether it differs from zero as derived from a two-tailed test. The results were obtained with phylogenetic linear regression by phylolm models on a single maximum clade credibility phylogenetic tree.Full size imageAllometry of appendagesAllen’s hypothesis3 implies that the length of animal’s appendages varies with temperature in relative (not absolute) terms, thus when asking how the appendage length vary across temperature gradient, we always need to control for body size. Phylogenetic log-log regression models revealed that body mass explains 72.7% and 72.5% of variance in beak and tarsus length (estimated with r-squared of models shown in Fig. 2a and Supplementary Fig. 3a), respectively, confirming that the evolution of absolute avian appendage size is substantially constrained by body size. These null allometric models predict that log-transformed beak length (Fig. 2a) and tarsus length (Fig. 3a) scale with log-transformed body mass in a linear manner:$${{{{rm{log }}}}}_{e}left({{{{rm{Beak}}}}},{{{{rm{Length}}}}}right)=1.4345+0.3362,{{{{rm{log }}}}}_{e}{{{{rm{Body}}}}},{{{{rm{Mass}}}}}$$
(1)
$${{{{rm{log }}}}}_{e}left({{{{rm{Tarsus}}}}},{{{{rm{Length}}}}}right)=2.1141+0.2883,{{{{rm{log }}}}}_{e}{Body},{{{{rm{Mass}}}}}$$
(2)
Fig. 2: Global test of Allen’s rule on avian beak length across 9962 (99.7%) bird species.The null allometric model (a) used to scale the absolute (log-transformed) beak length with log body size, the residuals from which were used as the relative beak length. Distribution of relative beak length across species geographic ranges (b). Model selection procedure for predicting log beak length (c), involving models with log body mass and either of six temperature measures within species geographic ranges included as fixed and interaction terms; AIC—Akaike Information Criterion, r2—coefficient of determination. An exemplar Allen’s model (d) showing increasing beak length with max temperature of all months, while controlling for body size as fixed term. An exemplar model with interaction of body size and max temperature of all months (e) illustrating how Allen’s rule operates across steeping quantiles of body size (left) and how allometry varies across quantiles of temperature (right). See Supplementary Fig. 2 for surrogate models based on the other temperature measures (evaluated in c). The p values refer to the significance of model’s fixed (d) or interaction terms (e) derived from two-tailed tests. The shaded area around the trend line is simple shading to facilitate reading. The results were obtained with phylogenetic linear regression by phylolm models on a single maximum clade credibility phylogenetic tree.Full size imageFig. 3: Global test of Allen’s rule on avian tarsus length across 9962 (99.7%) bird species.The null allometric model (a) used to scale the absolute (log-transformed) tarsus length with log body size, the residuals from which were used as the relative tarsus length. Distribution of relative tarsus length across species geographic ranges (b). Model selection procedure for predicting log tarsus length (c), involving models with log body mass and either of six temperature measures within species geographic ranges included as fixed and interaction terms; AIC—Akaike Information Criterion, r2—coefficient of determination. An exemplar Allen’s model (d) showing decreasing tarsus length with max temperature of all months, while controlling for body size as fixed term. An exemplar model with interaction of body size and max temperature of all months (e) illustrates how Allen’s rule operates across steeping quantiles of body size (left) and how allometry varies across steeping quantiles of temperature (right). See Supplementary Fig. 3 for surrogate models based on the other temperature measures (evaluated in c). The p values refer to the significance of model’s fixed (d) or interaction terms (e) derived from two-tailed tests. The shaded area around the trend line is simple shading to facilitate reading. The results were obtained with phylogenetic linear regression by phylolm models on a single maximum clade credibility phylogenetic tree.Full size imagethe normalized formulas of which give us the logarithmic equations:$${{{{rm{Beak}}}}},{{{{rm{Length}}}}}={4.1975,{{{{rm{Body}}}}},{{{{rm{Mass}}}}}}^{0.3362}$$
(3)
$${{{{rm{Tarsus}}}}},{{{{rm{Length}}}}}={8.2821,{{{{rm{Body}}}}},{{{{rm{Mass}}}}}}^{0.2883}$$
(4)
Because these allometric plots (Figs. 2a and 3a) relate the length of the appendage (one dimensional linear measure) to the body mass (three-dimensional volumetric measure) it means that the size of appendages would scale isometrically (proportionally) with the body size if the allometric coefficient was 0.3333. Thus, beak length equals to body mass to a power of 0.3362 means that the beak elongates almost exactly proportionally with body size. However, tarsus length equals to body mass to a power of 0.2883 means that the extent to which tarsus elongates with body size is slightly more pronounced in smaller species and weaker in larger species.These allometric relationships have implications for how we interpret subsequent patterns. For example, consider a species that experience a temporal increase in temperature, or invades a warmer climate. Then, if only Bergmann’s rule is operating (and in the absence of other confounds), a gradual decrease in body size should result in a proportional decrease in absolute beak length, and a gradually larger decrease in absolute tarsus length. Conversely, if species follow only Allen’s rule (and not follow Bergmann’s rule), then the increase in beak length should be similar between larger and smaller species, while the increase in tarsus length should be weaker in larger species and stronger in smaller species. Thus, Allen’s assumption that the increase in the ratio of body width to body length is steeper in larger species3, should not be a direct effect of the allometric rules, as appendages tend to increase proportionally with body size (beak) or increase milder at larger body sizes (tarsus).Allen’s ruleAfter excluding the effect of allometry, relative beak length is still tightly associated with phylogeny (Supplementary Table 1), while showing an impressive geographic variation (Fig. 2b). Our phylogenetic analysis concurs with an array of existing studies12,16,17,44,45 that found that the length of avian beak follows Allen’s rule, and is a general pattern across birds as a whole. Among our models predicting beak length, those with temperature variables among fixed terms are more informative than the null allometric model, where log body mass (allometry) is put as sole predictor (Fig. 2c). Most of the temperature variables also predict beak length better than latitude (Fig. 2c), again confirming the thermoregulatory basis of the observed pattern. Each of the temperature variables are positively associated with longer beaks (Fig. 2d and Supplementary Fig. 2a), which remains in agreement with Allen’s rule.Some studies have reported the ambiguous46 or very weak16 Allen’s pattern for avian legs. While relative tarsus length is also well conserved in avian phylogeny (Supplementary Table 1) and shows a high geographic variation (Fig. 3b), surprisingly, our global phylogenetic analysis indicates that avian tarsus length follows the inverse of Allen’s rule. Among models explaining tarsus length, those with temperature variables are better than the null allometric model (Fig. 3c). However, these models indicate a negative correlation—thereby shorter tarsi are associated with warmer temperatures (Fig. 3d).Allen’s vs Bergmann’s rule in allometryOur analyses support the hypothesis that the way in which avian appendages size varies across temperature regimes, depends on body size and vice versa. First, among models of beak length, those with an interaction of body size and the temperature consistently perform better than models without that interaction (Fig. 2c). The interaction of temperature and body size loads positively on beak length, indicating that larger-bodied species show stronger increases in beak size with temperature (Fig. 2e, left plot). Notably, beak length does not co-vary with temperature in the smallest birds (Fig. 2e, left plot), which is in agreement with Allen’s speculations3 that being smaller reduces the need to develop elongated appendages in hot climates, as effective heat exchange is already enabled through small body size (according to Bergmann’s rule). The positive interaction between body size and temperature also indicates that the higher the temperature, the steeper the allometric relationship between beak size and body size (Fig. 2e, right plot), meaning that in warmer climates beak size increases more strongly with body size than in colder climates, exactly as Allen hypothesized.An interaction between body size and temperature is also consistently supported in models of tarsus length (Fig. 3c). This interaction has strong positive effect on tarsus length, thereby reversing the trend by which tarsus shortens with temperature (Fig. 3e, left plot). This means that despite the overall decrease of tarsus size with temperature in smaller birds (the inverse of Allen’s rule), the opposite is true for larger birds that show increasing tarsus size with temperature (Fig. 3e, left plot). The interaction holds regardless of the temperature measure examined (Supplementary Fig. 3b, upper row), even if those previously did not co-vary with tarsus length when included as simple independent term with body size (Supplementary Fig. 3a). The case of larger birds thus fits Allen’s rule, and agrees with Allen’s further speculations3 that appendages are more likely to increase in larger- than in smaller-bodied animals. However, Allen did not predict the possibility of shortening appendages toward hot temperatures as seen in small birds. Given the extent of our sampling, the effect of shortening tarsi toward the equator in small-bodied species is presumably not an artefact, but relies on yet unknown mechanisms (possibly unrelated to thermoregulation). Nevertheless, if there is an evolutionary pressure to develop a smaller tarsus in hot climates, the increased thermoregulatory needs of larger-bodied species possibly overwhelm this selective process. This may be because large species acquire higher heat loads when the ambient temperature is hot, hence necessitating the development of longer legs as cooling organs. As with beak size, the interaction also indicates substantial changes in allometry, with much more millimeters of tarsus per each gram of body in warm conditions compared to cold conditions (Fig. 3e, right plot).Our analyses also support the mirror scenario, that the extent to which body size decreases with temperature (Bergmann’s rule) depends on the length of appendage. In models predicting body size, the temperature does not interact with relative beak length (Supplementary Fig. 4a), but interacts with tarsus length (Supplementary Fig. 4b). This interaction indicates that the strongest shrinkage in body size with temperature occurs in shorter-legged birds, while in longer-legged birds body size increases with temperature (inverse Bergmann’s rule). This again supports Allen’s speculations that variation in body shape allows birds to evolve body sizes less restricted (or even unrestricted) to environmental temperature. Thus, the results support the theory that Bergmann’s and Allen’s rules are two distinct, albeit analogous strategies to deal with thermoregulation.Allen’s vs Bergmann’s rule in climatic adaptationsOur analysis shows that the interactions of body size (Bergmann’s rule), beak length and tarsus length (Allen’s rule) predict the thermal environment across birds (e.g. the max temperatures of all months across species ranges, Fig. 4a). As with body size and shape, the temperatures experienced by species within their geographic ranges are finely conserved in the avian phylogeny (Supplementary Table 1), suggesting that thermal preferences of avian species have been established through evolutionary history. Evolution of these preferences then occurred when temperature changes affected their native environments (thus causing extinctions or adaptations), or when birds invaded novel environments (thus adapting to newly-encountered climates). Log-transformed body mass, relative beak length and tarsus length clearly predict the species ambient temperature (Fig. 4b), suggesting that the phenotype changes as animals adapt to suit different climates. However, of particular note is that the addition of an interaction between body size and relative beak length substantially improves model performance (Fig. 4b). This interaction shows that for longer-beaked birds, temperature associations are unrelated to body size, but the shorter the beak, the more pronounced is the shrinking in body size in warmer temperatures (Fig. 4c, left plot). In the case of smaller-bodied birds, the adaptation to different temperatures is independent of beak length, but with larger birds, the adaptation to warmer temperatures is more likely associated with elongated beaks (Fig. 4c, right plot). These results indicate that living in warmer temperatures tends to be associated either with smaller body size (Bergmann’s rule) or longer beak (Allen’s rule), rather than both rules simultaneously, thus again supporting the hypothesis of an evolutionary compromise between shifts in body size and shape as alternative adaptations to thermal environment.Fig. 4: Global test for avian adaptation to maximum temperature across all months by shifts in body size (Bergmann’s rule) and appendage size (Allen’s rule) across 9962 (99.7%) avian species.Distribution of environmental temperature across species geographic ranges (a). Model selection procedure for predicting max temperature all months (b), involving models with different combinations of log body mass, relative beak and tarsus length as fixed and interaction terms; AIC—Akaike Information Criterion, r2—coefficient of determination. Exemplar models with two-way interaction of body size and relative beak length (c) or tarsus length (d) illustrate how Bergmann’s rule operate across steeping quantiles of relative appendage length (left plots) and how Allen’s rule operate across steeping quantiles of body size (right plots). An exemplar model with two-way interaction of relative beak length and tarsus length (e) illustrates how Allen’s rule based on the relative length of one appendage operates across steeping quantiles of the relative length of second appendage. An exemplar model with three-way interaction of log body mass, relative beak and tarsus length (f) illustrates how shifts in body size and two measures of body shape depend on each other when animals adapt to novel climates; the trend lines indicate relationships between y and x1 (axes) across combinations of min and max values of x2 and x3 (colors); see also Supplementary Fig. 5 for more detailed visualization of the model f. The p values refer to significance of two-way (c–e) and three-way (f) interaction terms derived from two-tailed tests. The shaded area around the trend line is simple shading to facilitate reading. Obtained with phylogenetic linear regression by phylolm models on a single maximum clade credibility phylogenetic tree.Full size imageThe interaction of body size and relative tarsus length also substantially improves the model predicting ambient temperature of the species (Fig. 4b). This interaction indicates that living in warmer climates is associated with smaller body size (Bergmann’s rule) only in shorter-legged birds, while in longer-legged birds the environmental temperature increases with body size (inverse Bergmann’s rule; Fig. 4d, left plot). Simultaneously, the avian environmental temperature increases with tarsus length (Allen’s rule) only in larger species, while the opposite is true for smaller species (Fig. 4d, right plot). This suggests that larger-legged avian lineages may be resistant to Bergmann’s rule and become larger when habituating to warm climates, while shorter- and average-legged birds become smaller with temperature, as predicted by Bergmann’s rule. These findings again converge with Allen’s speculations on trade-off in the evolution of body size and appendage length in relation to temperature.We found that the length of the two different appendages—beak and tarsus—show independent evolutionary patterns (Fig. 4e). The environmental temperature of a species increases with beak length independently from tarsus length, and decreases with tarsus length independently from beak length (Fig. 4e). These outcomes reject the possibility of an evolutionary compromise in climatic adaptation of two types of appendages, at least when we do not control for body size (Bergmann’s rule) as additional type of climatic adaptation.Finally, the model with a three-way interaction between body size, relative beak and tarsus length predicting temperature performs the best among all considered candidate models (Fig. 4b) and this interaction is statistically significant (Fig. 4f), suggesting that evolutionary adaptation to novel climates depends on various configurations of body size, beak, and tarsus length. This model indicates various Bergmann’s rule slopes across different settings of body shape (Fig. 4f, top-left). Namely, the steepest decrease in environmental temperature with body size (i.e. strongest Bergmann’s rule) is observed in smaller-billed and smaller-legged birds (Fig. 4f, top-left, brown trend line), whereas in longer-billed and shorter-legged birds (Fig. 4f, top-left, green trend line) body size is not associated with environmental temperature. This model also indicates that in shorter-billed, longer-legged birds (Fig. 4f, top-left, purple trend line) body size increases across temperature gradient (inverse Bergmann’s pattern). This thus strengthens the support for Allen’s theory that having bodies with elongated appendages may enable species to circumvent or even reverse Bergmann’s pattern; whereas compact bodies are more prone to decrease in size with temperature in order to deal with overheating in warm climates. Counteracting this argument, however, is that longer-billed and longer-legged birds show (moderate) typical Bergmann’s pattern (Fig. 4f, leftmost plot, bluish trend line).The three-way interaction model also shows other mixtures of expected and unexpected results. For example, the strongest increase in environmental temperature with beak length occurs in larger-bodied and shorter-legged birds (Fig. 4f, top-right plot, orange trend line), which clearly suggests a trade-off in evolution of body size and beak length and a similar trade-off in the evolution of the two types of appendages, presumably reflecting different adaptive responses for thermoregulation. However, a similar increase in beak length also occurs in tiny-bodied and longer-legged birds (Fig. 4f, top-right plot, blue trend line), which stands in contrast to this trade-off hypothesis. Likewise, the steepest increase in environmental temperature with tarsus length (Allen’s rule) occurs in larger-bodied and shorter-billed birds (Fig. 4f, bottom plot, pink trend line), again suggesting a compromise scenario, with elongated tarsus evolving as thermoregulatory organ to compensate for insufficient heat exchange due to large body and small beak. It also suggests that, in large birds, having a short beak in hot climates requires longer tarsi (Fig. 4f, top-right, pinkish and orange trend lines) and vice versa (Fig. 4f, bottom plot, rose and yellowish trend lines), indicating that in large species, the summarized length of two types of appendages is important for thermoregulation. However, by contrast, it seems that in small bodied species, beak and tarsi length evolved in a correlated way (Fig. 4f, top-right and bottom plots, green and blue trend lines) across environmental temperature (occurrences in warmer temperatures are associated with simultaneously both longer beaks and tarsi, or else simultaneously shorter beaks and tarsi). This may indicate a general tendency to correlated evolution of relative beak and tarsus lengths, perhaps for functional reasons, e.g. longer beaks may allow long-legged birds to explore substrate more efficiently, as longer necks also do47.Allen’s vs Bergmann’s rules in causal modelsOur hypothesis consequently holds within phylogenetic path analysis, where the best causal models integrate Bergmann’s and Allen’s rules to explain both the size of avian appendages (Fig. 5a) and the avian thermal environment (here, maximum temperature across all months) (Fig. 5b). The best model predicting beak and tarsus length includes the causal effect of temperature on body size (Bergmann’s rule) and then body size on beak and tarsus length (allometry), as well as the direct effect of temperature on the size of appendages (Allen’s rule) (Fig. 5a). This joint Bergmann’s and Allen’s model is substantially better than the model assuming that temperature does not affect body size before scaling for the length of appendages (Fig. 5a, Allen’s rule only). The combined Bergmann’s and Allen’s model is also better than one assuming no direct effect of temperature on appendages (Fig. 5a, Bergmann’s rule only). This again indicates that how the length of avian appendages co-varies with the ambient temperature partially depends on how avian body size co-varies with temperature, yielding results aligned with the trade-off hypothesis. This notably argues against the possibility that the increase in the length of appendages (relative to body size) with temperature is an artefact of decreased body sizes at hot temperatures (see26). However, interestingly, the model including only Allen’s rule (and allometry) explains the length of appendages with similar accuracy to the model with only Bergmann’s rule (Fig. 5a).Fig. 5: Phylogenetic path analysis with responses of the length of avian appendages (beak and tarsus) (a) and the maximum temperature of all months within species range (b) across 9962 species (99.7% of global community).In both cases model candidates include different combinations of allometry (relationship between the length of appendages and body size), Bergmann’s rule (the relationship between body size and the temperature) and Allen’s rule (relationship between the length of appendages and temperature). ∆CIC—delta C statistic Information Criterion. The results were obtained with phylogenetic path analysis by using phylopath models on a single maximum clade credibility phylogenetic tree and scaled covariates (mean = 0 ± 1 SD) to compare their effect sizes (see numbers on path diagrams).Full size imageThe best model predicting the temperature associations includes the indirect effect of body size on the length of appendages (allometry), and then the length of appendages on temperature (Allen’s rule), as well as the direct effect of body size on temperature (Bergmann’s rule) (Fig. 5b). These results again demonstrate that how the temperature varies across species ranges depends on both the size of body and appendages, suggesting that Bergmann’s and Allen’s rules describe two distinct evolutionary ways to cope with thermoregulation. Moreover, the similar performance of Allen’s model compared to Bergmann’s model (Fig. 5b) again suggests that shifts in the animal’s body size and shape represent roughly equally influential in the evolution of adaptations to novel climates.Excluding possible confounding factorsTo ensure the reliability of our findings, first we show that when explaining the phenotype (Supplementary Figs. 2b and 3b), or the temperature within species geographic ranges (Supplementary Fig. 6), the main results remain consistent whichever of the five temperature measures is included. Second, despite the fact that the relationships with relative length of appendages and the experienced temperatures are strongest in resident birds, followed by partial- and full- migrants, our results still hold when accounting for these three categories of avian migratory habits (Supplementary Fig. 7); and the compromise scenario remains similar in each of these groups independently (Supplementary Figs. 8–10). It aligns with previous studies22,46, which found that ecogeographical rules are valid regardless of variation in avian migratory habits. However, it is worth to notice that the most prominent trade-offs are found in resident species (in case of explaining environmental temperature, see Supplementary Fig. 10) or in partial migrants (in case of explaining beak length, see Supplementary Fig. 8). Third, the trade-offs in thermoregulatory strategies also hold after controlling for geographic range size (Supplementary Fig. 11) and remains quantitatively (Supplementary Figs. 12–13) or qualitatively (Supplementary Fig. 14) stable across the gradient of endemic-cosmopolitan species. Thus, even if ecogeographical rules operate within widespread species (across distanced populations, as well documented9,10,11,12,13), this does not appear to influences the results of our cross-species analysis. Fourth, the predictions of temperature within species geographic ranges are also not specific to the way by which we account for the allometry of appendages (by using residual appendage length). Parallel analyses with ratios of appendage length to body mass (Supplementary Fig. 15) or principal components of all phenotypic traits (Supplementary Fig. 16) give qualitatively similar outcomes. Fifth, we also show that results of both phylogenetic regression (Supplementary Fig. 17) and phylogenetic path models (Supplementary Fig. 18) remain consistently valid across 100 randomly chosen phylogenetic trees32, mitigating concerns regarding phylogenetic uncertainty influencing our results.Notably, there is a wider list of important ecological factors constraining or favoring variation in body size and shape, e.g. tropic levels or foraging techniques21,23,47,48, although they are also themselves constrained by phylogeny to some extent, which we control for. Nevertheless, we believe that it is likely that these constraints influenced (or were influenced by) the Bergmann-Allen trade-off. Understanding of this issue would benefit from a deeper dive into the relationship between climatic, phenotypic and ecological variation across animals.Our findings in the context of eco-evolutionary processes driven by climateTo the best of our knowledge, this study is the largest (taxonomically and geographically) simultaneous test of ecogeographical rules and it provides a first empirical evidence for a trade-off in the evolution of body size (Bergmann’s rule2) and the size of appendages (Allen’s rule3) across global temperature gradients. Our results confirm what Allen3 speculated—the larger the body, the stronger the increase in appendage size with temperature; and the larger the appendages, the milder the decrease in body size with temperature. Thus, the evolution of body size under temperature regimes likely depends on the size of appendages and, on the other hand, the extent to which temperature drives the size of appendages depends on body size. This means that these two thermoregulatory adaptations are not independent of each other, but the phenotype has at least two ways to adapt to novel climates, i.e. by the shifts in body size or the shifts in the size of appendages (or both to a lesser extent).The evolution of appendages (e.g. avian beaks49) was a dynamic process believed to overtake the changes in body size across evolutionary time50. Our analyses do not indicate, however, that shifts in body size have been more frequent than shifts in appendage size (or vice versa), at least not because of thermoregulation. Rather, they indicate that shifts in body size and shape are intertwined through avian evolutionary history, agreeing with the theory that animals select the most convenient strategy of thermoregulation to maintain functional traits of its phenotype. For functional reasons animal lineages tend to increase in body size over evolutionary time (Cope’s rule43), thus it is not surprising that strategies allowing species to maintain/develop larger bodies (i.e. over-increase in appendage size) are to be expected evolutionarily. On the other hand, some lineages may be constrained in appendage size (e.g. to forage21,47 or communicate23 effectively), hence those may favor the shifts in body size to reconcile optimal thermoregulation with a desired functionality.We found that the compromise in thermoregulatory strategies may also involve two distinct types of appendages, here beak versus tarsus. However, this is true only for larger-bodied species (see Fig. 4f top-right and bottom plots, trends for large bodies), that are more likely to acquire higher heat loads in warmer environments, thus the summarized size of many appendages may be for them crucial to disperse heat loads. Both beak and legs have been confirmed to act as key regions of heat transfer on the avian body37,38,51,52, thus both may be sensitive to thermal conditions when body size is too large to deal alone with too hot temperatures. Yet, in small-bodied species both appendages seem to evolve in concerted way across temperature gradients, and this may be in a way that conforms with Allen’s rule or not (see Fig. 4f, top-right and bottom plots, trends for small bodies), indicating that the small body ensures good temperature exchange in hot climates, thus the evolution of appendages in these species may be correlated, but independent of thermoregulatory selection pressures.It is worthwhile emphasizing that apart from shifts in body size and shape, many other elements combine to help birds meet their thermoregulatory requirements53, e.g. through variation in insulation (feathers)54, coloration55,56 metabolism57, blood circulation58 or behavior59,60,61. Extrapolating our results, these thermoregulatory strategies might also co-evolve under a trade-off to ensure optimal thermoregulation along with desired functionality. This is presumably a reason for the relatively low performance of our models; e.g. physical phenotype explains up to only 20% of the variance in ambient temperature (Fig. 4b, upper model), therefore unexplained variance must be attributed to other thermoregulatory strategies.In this study, we demonstrate that Allen’s rule may be attributed to the varying allometric functions across temperature gradients. Although logical and argued elsewhere26, it has never been addressed by any empirical research. Our findings clearly indicate the importance of considering body mass as both a fixed and interaction term in studies of Allen’s rule, but also might suggest that ambient temperature should be included in other allometric studies of animals’ morphology. That said, temperature explains very little of the variance in the size of appendages compared to body size (Figs. 2c and 3c), thus thermal conditions are unlikely to be a very crucial confounding factor for allometry in comparative analyses.In this study, we empirically confirm for the first time an evolutionary compromise theory that was first proposed almost 150 years ago3– the evolution of body size and appendages are two distinct and interacting ways to cope with thermoregulation. This may explain why many studies fail to detect Allen’s or Bergmann’s rules independently which has led to questioning of the generality of these ecogeographical patterns13,24,25. Here, our findings suggest that Bergmann’s and Allen’s rules should not necessarily be considered in isolation. We believe that these thermoregulatory strategies might intertwine through the evolutionary history of animals, as the evolution of phenotype possibly interacts to confound ecogeographical rules to evolve functional traits. This explanation also highlights the diverse mechanisms that animals may employ to expand across the world’s multiple environments. It also raises the speculation that with observed and future anticipated warming of Earth’s climate, we should expect mainly large animals to elongate in appendages, while mainly compact-bodied animals to shrink in size. More

Tilman, D., Cassman, K. G., Matson, P. A., Naylor, R. & Polasky, S. Agricultural sustainability and intensive production practices. Nature 418, 671–677 (2002).Article 
CAS 
PubMed 

Google Scholar 
Tilman, D. et al. Forecasting agriculturally driven global environmental change. Science 292, 281–284 (2001).Article 
CAS 
PubMed 

Google Scholar 
IPBES: Summary for Policymakers. In The Assessment Report on Pollinators, Pollination and Food Production (eds Potts, S. G. et al.) (IPBES, 2016).Potts, S. G. et al. Safeguarding pollinators and their values to human well-being. Nature 540, 220–229 (2016).Article 
CAS 
PubMed 

Google Scholar 
Sgolastra, F. et al. Synergistic mortality between a neonicotinoid insecticide and an ergosterol-biosynthesis-inhibiting fungicide in three bee species. Pest Manag Sci. 73, 1236–1243 (2016).Article 
PubMed 

Google Scholar 
Whitehorn, P. R., O’Connor, S., Wackers, F. L. & Goulson, D. Neonicotinoid pesticide reduces bumble bee colony growth and queen production. Science 336, 351–352 (2012).Article 
CAS 
PubMed 

Google Scholar 
Rundlöf, M. et al. Seed coating with a neonicotinoid insecticide negatively affects wild bees. Nature 521, 77–80 (2015).Article 
PubMed 

Google Scholar 
Woodcock, B. et al. Impacts of neonicotinoid use on long-term population changes in wild bees in England. Nat. Commun. 7, 12459 (2016).Stuligross, C. & Williams, N. Past insecticide exposure reduces bee reproduction and population growth rate. Proc. Natl Acad. Sci. USA 118, e2109909118 (2021).Stanley, D. A. et al. Neonicotinoid pesticide exposure impairs crop pollination services provided by bumblebees. Nature 528, 548–550 (2015).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Tamburini, G. et al. Fungicide and insecticide exposure adversely impacts bumblebees and pollination services under semi-field conditions. Environ. Int. 157, 106813 (2021).Article 
CAS 
PubMed 

Google Scholar 
Sponsler, D. B. et al. Pesticides and pollinators: a socioecological synthesis. Sci. Total Environ. 662, 1012–1027 (2019).Article 
CAS 
PubMed 

Google Scholar 
Meehan, T. D., Werling, B. P., Landis, D. A. & Gratton, C. Agricultural landscape simplification and insecticide use in the Midwestern United States. Proc. Natl Acad. Sci. USA 108, 11500–11505 (2011).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Nicholson, C. C. & Williams, N. M. Cropland heterogeneity drives frequency and intensity of pesticide use. Environ. Res. 16, 074008 (2021).CAS 

Google Scholar 
Böhme, F., Bischoff, G., Zebitz, C. P. W., Rosenkranz, P. & Wallner, K. Pesticide residue survey of pollen loads collected by honeybees (Apis mellifera) in daily intervals at three agricultural sites in South Germany. PLoS ONE 13, e0199995 (2018).Larsen, A. E. & Noack, F. Impact of local and landscape complexity on the stability of field-level pest control. Nat. Sustain. 4, 120–128 (2021).Article 

Google Scholar 
Botías, C. et al. Neonicotinoid residues in wildflowers, a potential route of chronic exposure for bees. Environ. Sci. Technol. 49, 12731–12740 (2015).Article 
PubMed 

Google Scholar 
Krupke, C. H., Holland, J. D., Long, E. Y. & Eitzer, B. D. Planting of neonicotinoid-treated maize poses risks for honey bees and other non-target organisms over a wide area without consistent crop yield benefit. J. Appl. Ecol. 54, 1449–1458 (2017).Article 
CAS 

Google Scholar 
Wintermantel, D. et al. Neonicotinoid-induced mortality risk for bees foraging on oilseed rape nectar persists despite EU moratorium. Sci. Total Environ. 704, 135400 (2020).Article 
CAS 
PubMed 

Google Scholar 
Krupke, C. H., Hunt, G. J., Eitzer, B. D., Andino, G. & Given, K. Multiple routes of pesticide exposure for honey bees living near agricultural fields. PLoS ONE 7, e29268 (2012).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Long, E. Y. & Krupke, C. H. Non-cultivated plants present a season-long route of pesticide exposure for honey bees. Nat. Commun. 7, 11629 (2016).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
David, A. et al. Widespread contamination of wildflower and bee-collected pollen with complex mixtures of neonicotinoids and fungicides commonly applied to crops. Environ. Int. 88, 169–178 (2016).Article 
CAS 
PubMed 

Google Scholar 
Heinrich, B. The foraging specializations of individual bumblebees. Ecol. Monogr. 46, 105–128 (1976).Article 

Google Scholar 
Bolin, A., Smith, H. G., Lonsdorf, E. V. & Olsson, O. Scale-dependent foraging tradeoff allows competitive coexistence. Oikos 127, 1575–1585 (2018).Article 

Google Scholar 
Cresswell, J. E., Osborne, J. L. & Goulson, D. An economic model of the limits to foraging range in central place foragers with numerical solutions for bumblebees. Ecol. Entomol. 25, 249–255 (2000).Article 

Google Scholar 
Rundlöf, M. et al. Flower plantings support wild bee reproduction and may also mitigate pesticide exposure effects. J. Appl. Ecol. 59, 2117–2127 (2022).Article 

Google Scholar 
Graham, K. K. et al. Identities, concentrations, and sources of pesticide exposure in pollen collected by managed bees during blueberry pollination. Sci. Rep. 11, 16857 (2021).Centrella, M. et al. Diet diversity and pesticide risk mediate the negative effects of land use change on solitary bee offspring production. J. Appl. Ecol. 57, 1031–1042 (2020).Article 
CAS 

Google Scholar 
De Palma, A. et al. Ecological traits affect the sensitivity of bees to land-use pressures in European agricultural landscapes. J. Appl. Ecol. 52, 1567–1577 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
Sponsler, D. B. & Johnson, R. M. Mechanistic modeling of pesticide exposure: the missing keystone of honey bee toxicology. Environ. Toxicol. Chem. 36, 871–881 (2017).Article 
CAS 
PubMed 

Google Scholar 
Holzschuh, A., Dormann, C. F., Tscharntke, T. & Steffan-Dewenter, I. Mass-flowering crops enhance wild bee abundance. Oecologia 172, 477–484 (2013).Article 
PubMed 

Google Scholar 
McArt, S. H., Fersch, A. A., Milano, N. J., Truitt, L. L. & Böröczky, K. High pesticide risk to honey bees despite low focal crop pollen collection during pollination of a mass blooming crop. Sci. Rep. 7, 46554 (2017).Sanchez-Bayo, F. & Goka, K. Pesticide residues and bees—a risk assessment. PLoS ONE 9, e94482 (2014).Article 
PubMed 
PubMed Central 

Google Scholar 
Zioga, E., Kelly, R., White, B. & Stout, J. C. Plant protection product residues in plant pollen and nectar: a review of current knowledge. Environ. Res. 189, 109873 (2020).Article 
CAS 
PubMed 

Google Scholar 
The European Green Deal (European Commission, 2019).More, S. J., Auteri, D., Rortais, A. & Pagani, S. EFSA is working to protect bees and shape the future of environmental risk assessment. EFSA J. 19, e190101 (2021).Schmolke, A. et al. Assessment of the vulnerability to pesticide exposures across bee species. Environ. Toxicol. Chem. 40, 2640–2651 (2021).Article 
CAS 
PubMed 

Google Scholar 
Rollin, O. et al. Differences of floral resource use between honey bees and wild bees in an intensive farming system. Agric. Ecosyst. Environ. 179, 78–86 (2013).Article 

Google Scholar 
Persson, A. S. & Smith, H. G. Seasonal persistence of bumblebee populations is affected by landscape context. Agric. Ecosyst. Environ. 165, 201–209 (2013).Article 

Google Scholar 
Samuelson, A. E., Schürch, R. & Leadbeater, E. Dancing bees evaluate central urban forage resources as superior to agricultural land. J. Appl. Ecol. 59, 79–88 (2022).Article 

Google Scholar 
Milner, A. M. & Boyd, I. L. Toward pesticidovigilance. Science 357, 1232–1234 https://doi.org/10.1126/science.aan2683 (2017).Nowell, L. H., Norman, J. E., Moran, P. W., Martin, J. D. & Stone, W. W. Pesticide toxicity index—a tool for assessing potential toxicity of pesticide mixtures to freshwater aquatic organisms. Sci. Total Environ. 476–477, 144–157 (2014).Article 
PubMed 

Google Scholar 
Mullin, C. A., Frazier, M., Frazier, J. L., Ashcraft, S. & Simonds, R. High levels of miticides and agrochemicals in North American apiaries: implications for honey bee health. PLoS ONE 5, 9754 (2010).Article 

Google Scholar 
Pettis, J. S. et al. Crop pollination exposes honey bees to pesticides which alters their susceptibility to the gut pathogen Nosema ceranae. PLoS ONE 8, e70182 (2013).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Végh, R., Sörös, C., Majercsik, N. & Sipos, L. Determination of pesticides in bee pollen: validation of a multiresidue high-performance liquid chromatography-mass spectrometry/mass spectrometry method and testing pollen samples of selected botanical origin. J. Agric. Food Chem. 70, 1507–1515 (2022).Article 
PubMed 

Google Scholar 
Park, M. G., Blitzer, E. J., Gibbs, J., Losey, J. E. & Danforth, B. N. Negative effects of pesticides on wild bee communities can be buffered by landscape context. Proc. R. Soc. B 282, 20150299 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
Graham, K. K. et al. Pesticide risk to managed bees during blueberry pollination is primarily driven by off-farm exposures. Sci. Rep. 12, 7189 (2022).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Yourstone, J., Karlsson, M., Klatt, B. K., Olsson, O. & Smith, H. G. Effects of crop and non-crop resources and competition: high importance of trees and oilseed rape for solitary bee reproduction. Biol. Conserv. 261, 109249 (2021).Persson, A. S., Mazier, F. & Smith, H. G. When beggars are choosers—how nesting of a solitary bee is affected by temporal dynamics of pollen plants in the landscape. Ecol. Evol. 8, 5777–5791 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Wood, T. J., Holland, J. M. & Goulson, D. Providing foraging resources for solitary bees on farmland: current schemes for pollinators benefit a limited suite of species. J. Appl. Ecol. 54, 323–333 (2016).Garthwaite, D. et al. Collection of Pesticide Application Data in View of Performing Environmental Risk Assessments for Pesticides (EFSA, 2017).de Oliveira, R. C., Nascimento Queiroz, S. C., Pinto da Luz, C. F., Silveira Porto, R. & Rath, S. Bee pollen as a bioindicator of environmental pesticide contamination. Chemosphere 163, 525–534 (2016).Article 
PubMed 

Google Scholar 
Arena, M. & Sgolastra, F. A meta-analysis comparing the sensitivity of bees to pesticides. Ecotoxicology 23, 324–334 (2014).Article 
CAS 
PubMed 

Google Scholar 
Douglas, M. R., Sponsler, D. B., Lonsdorf, E. V. & Grozinger, C. M. County-level analysis reveals a rapidly shifting landscape of insecticide hazard to honey bees (Apis mellifera) on US farmland. Sci. Rep. 10, 797 (2020).Commission Implementing Regulation (EU) 2021/2081 of 26 November 2021 concerning the non-renewal of approval of the active substance indoxacarb, in accordance with Regulation (EC) No 1107/2009 of the European Parliament and of the Council concerning the placing of plant protection products on the market, and amending Commission Implementing Regulation (EU) No 540/2011 (EUR-Lex, 2021); http://data.europa.eu/eli/reg_impl/2021/2081/ojCommission Implementing Regulation (EU) 2020/23 of 13 January 2020 concerning the non-renewal of the approval of the active substance thiacloprid, in accordance with Regulation (EC) No. 1107/2009 of the European Parliament and of the Council concerning the placing of plant protection products on the market, and amending the Annex to Commission Implementing Regulation (EU) No 540/2011 (EUR-Lex, 2020); http://data.europa.eu/eli/reg_impl/2020/23/ojCommission Implementing Regulation (EU) 2018/783 of 29 May 2018 amending Implementing Regulation (EU) No 540/2011 as regards the conditions of approval of the active substance imidacloprid (EUR-Lex, 2018); http://data.europa.eu/eli/reg_impl/2018/783/ojHerbertsson, L., Jonsson, O., Kreuger, J., Smith, H. G. & Rundlöf, M. Scientific note: imidacloprid found in wild plants downstream permanent greenhouses in Sweden. Apidologie 52, 946–949 (2021).Article 

Google Scholar 
Tosi, S. et al. Long-term field-realistic exposure to a next-generation pesticide, flupyradifurone, impairs honey bee behaviour and survival. Commun. Biol. 4, 805 (2021).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Siviter, H. & Muth, F. Do novel insecticides pose a threat to beneficial insects?: novel insecticides harm insects. Proc. R. Soc. B 287, 20201265 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
EFSA. Guidance on the risk assessment of plant protection products on bees (Apis mellifera, Bombus spp. and solitary bees). EFSA J. 11, 3295 (2013).Guidance for Assessing Pesticide Risks to Bees (US EPA, 2014).Boyle, N. K. et al. Workshop on pesticide exposure assessment paradigm for non-apis bees: foundation and summaries. Environ. Entomol. 48, 4–11 (2019).Article 
PubMed 

Google Scholar 
EFSA. Analysis of the evidence to support the definition of specific protection goals for bumble bees and solitary bees. EFSA J. 19, EN-7125 (2022).Garibaldi, L. A. et al. Wild pollinators enhance fruit set of crops regardless of honey bee abundance. Science 339, 1608–1611 (2013).Article 
CAS 
PubMed 

Google Scholar 
Tscharntke, T., Grass, I., Wanger, T. C. & Westphal, C. Restoring biodiversity needs more than reducing pesticides. Trends Ecol. Evol. 37, 115–116 (2022).Article 
PubMed 

Google Scholar 
Topping, C. J. et al. Holistic environmental risk assessment for bees. Science 37, 897 (2021).Article 

Google Scholar 
Tsvetkov, N. et al. Chronic exposure to neonicotinoids reduces honey bee health near corn crops. Science 356, 1395–1397 (2017).Article 
CAS 
PubMed 

Google Scholar 
Jonsson, O., Fries, I. & Kreuger, J. Utveckling av Analysmetoder och Screening av Växtskyddsmedel i bin och Pollen (CKB, 2013).Sawyer, R. Pollen Identification for Beekeepers (Univ. Cardiff Press, 1981).IUPAC Pesticide Properties Data Base (Univ. of Hertfordshire, 2022).EFSA Scientific Committee & More, S.J. et al. Guidance on harmonised methodologies for human health, animal health and ecological risk assessment of combined exposure to multiple chemicals. EFSA J. 17, e05634 (2019).Martin, O. et al. Ten years of research on synergisms and antagonisms in chemical mixtures: a systematic review and quantitative reappraisal of mixture studies. Environ. Int. 146, 106206 (2021).Article 
CAS 
PubMed 

Google Scholar 
DiBartolomeis, M., Kegley, S., Mineau, P., Radford, R. & Klein, K. An assessment of acute insecticide toxicity loading (AITL) of chemical pesticides used on agricultural land in the United States. PLoS ONE 14, e0220029 (2019).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Test No. 213: Honeybees, Acute Oral Toxicity Test (OECD, 1998); https://doi.org/10.1787/9789264070165-enPrice, P. S. & Han, X. Maximum cumulative ratio (MCR) as a tool for assessing the value of performing a cumulative risk assessment. Int. J. Environ. Res. Public Health 8, 2212–2225 (2011).Article 
PubMed 
PubMed Central 

Google Scholar 
Bates, D., Mächler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 1–48 (2015).Oksanen, J. et al. vegan community ecology package version 2.6-2 (2022).Lenth, R. emmeans: Estimated marginal means, aka least-squares means (2022).Lüdecke, D., Ben-shachar, M. S., Patil, I. & Makowski, D. performance: an R package for assessment, comparison and testing of statistical models statement of need. J. Open Source Softw. 6, 3139 (2021).Nakagawa, S. & Schielzeth, H. A general and simple method for obtaining R2 from generalized linear mixed-effects models. Methods Ecol. Evol. 4, 133–142 (2013).Article 

Google Scholar 
Kendall, L. K. et al. The potential and realized foraging movements of bees are differentially determined by body size and sociality. Ecology 103, e3809 (2022).Parreño, M. A. et al. Critical links between biodiversity and health in wild bee conservation. Trends Ecol. Evol. 37, 309–321 (2022).Article 
PubMed 

Google Scholar  More

Semin, G. R., Scandurra, A., Baragli, P., Lanatà, A. & D’Aniello, B. Inter- and intra-species communication of emotion: Chemosignals as the neglected medium. Animals 9, 887 (2019).Article 
PubMed 
PubMed Central 

Google Scholar 
Désiré, L., Boissy, A. & Veissier, I. Emotions in farm animals: A new approach to animal welfare in applied ethology. Behav. Processes 60, 165–180 (2002).Article 
PubMed 

Google Scholar 
Briefer, E. F. & Le Comber, S. Vocal expression of emotions in mammals: mechanisms of production and evidence. J. Zool. 288, 1–20 (2012).Article 

Google Scholar 
Jardat, P. & Lansade, L. Cognition and the human–animal relationship: a review of the sociocognitive skills of domestic mammals toward humans. Anim. Cogn. 25, 369–384 (2022).Article 
PubMed 

Google Scholar 
Galvan, M. & Vonk, J. Man’s other best friend: domestic cats (F. silvestris catus) and their discrimination of human emotion cues. Anim. Cogn. 19, 193–205 (2016).Nawroth, C., Albuquerque, N., Savalli, C., Single, M. S. & McElligott, A. G. Goats prefer positive human emotional facial expressions. R. Soc. Open Sci. 5, 180491 (2018).Article 
ADS 
PubMed 
PubMed Central 

Google Scholar 
Proops, L., Grounds, K., Smith, A. V. & McComb, K. Animals remember previous facial expressions that specific humans have exhibited. Curr. Biol. 28, 1428-1432.e4 (2018).Article 
CAS 
PubMed 

Google Scholar 
Smith, A. V., Proops, L., Grounds, K., Wathan, J. & McComb, K. Functionally relevant responses to human facial expressions of emotion in the domestic horse (Equus caballus). Biol. Lett. 12, 20150907 (2016).Article 
PubMed 
PubMed Central 

Google Scholar 
Siniscalchi, M., D’Ingeo, S., Fornelli, S. & Quaranta, A. Lateralized behavior and cardiac activity of dogs in response to human emotional vocalizations. Sci. Rep. 8, 77 (2018).Article 
ADS 
PubMed 
PubMed Central 

Google Scholar 
Siniscalchi, M., D’Ingeo, S. & Quaranta, A. Orienting asymmetries and physiological reactivity in dogs’ response to human emotional faces. Learn. Behav. 46, 574–585 (2018).Article 
PubMed 

Google Scholar 
Smith, A. V. et al. Domestic horses (Equus caballus) discriminate between negative and positive human nonverbal vocalisations. Sci. Rep. 8, 13052 (2018).Article 
ADS 
PubMed 
PubMed Central 

Google Scholar 
Trösch, M. et al. Horses categorize human emotions cross-modally based on facial expression and non-verbal vocalizations. Animals 9, 862 (2019).Article 
PubMed 
PubMed Central 

Google Scholar 
Quaranta, A., D’ingeo, S., Amoruso, R. & Siniscalchi, M. Emotion recognition in cats. Animals 10, 1107 (2020).Nakamura, K., Takimoto-Inose, A. & Hasegawa, T. Cross-modal perception of human emotion in domestic horses (Equus caballus). Sci. Rep. 8, 8660 (2018).Article 
ADS 
PubMed 
PubMed Central 

Google Scholar 
Albuquerque, N. et al. Dogs recognize dog and human emotions. Biol. Lett. 12, 20150883 (2016).Article 
PubMed 
PubMed Central 

Google Scholar 
Briefer, E. F. Vocal contagion of emotions in non-human animals. Proc. R. Soc. B Biol. Sci. 285 (2018).Sabiniewicz, A., Tarnowska, K., Świątek, R., Sorokowski, P. & Laska, M. Olfactory-based interspecific recognition of human emotions: Horses (Equus ferus caballus) can recognize fear and happiness body odour from humans (Homo sapiens). Appl. Anim. Behav. Sci. 230, 105072 (2020).Article 

Google Scholar 
Baba, C., Kawai, M. & Takimoto-Inose, A. Are horses (Equus caballus) sensitive to human emotional cues?. Animals 9, 630 (2019).Article 
PubMed 
PubMed Central 

Google Scholar 
Brennan, P. A. & Kendrick, K. M. Mammalian social odours: attraction and individual recognition. Philos. Trans. R. Soc. B Biol. Sci. 361, 2061–2078 (2006).Saslow, C. A. Understanding the perceptual world of horses. Appl. Anim. Behav. Sci. 78, 209–224 (2002).Article 

Google Scholar 
Péron, F., Ward, R. & Burman, O. Horses (Equus caballus) discriminate body odour cues from conspecifics. Anim. Cogn. 17, 1007–1011 (2014).Article 
PubMed 

Google Scholar 
Krueger, K. & Flauger, B. Olfactory recognition of individual competitors by means of faeces in horse (Equus caballus). Anim. Cogn. 14, 245–257 (2011).Article 
PubMed 

Google Scholar 
Boissy, A., Terlouw, C. & Le Neindre, P. Presence of cues from stressed conspecifics increases reactivity to aversive events in cattle: evidence for the existence of alarm substances in urine. Physiol. Behav. 63, 489–495 (1998).Article 
CAS 
PubMed 

Google Scholar 
Vieuille-Thomas, C. & Signoret, J. P. Pheromonal transmission of an aversive experience in domestic pig. J. Chem. Ecol. 18, 1551–1557 (1992).Article 
CAS 
PubMed 

Google Scholar 
Siniscalchi, M., D’Ingeo, S. & Quaranta, A. The dog nose ‘KNOWS’ fear: Asymmetric nostril use during sniffing at canine and human emotional stimuli. Behav. Brain Res. 304, 34–41 (2016).Article 
PubMed 

Google Scholar 
Calvi, E. et al. The scent of emotions: A systematic review of human intra- and interspecific chemical communication of emotions. Brain Behav. 10 (2020).de Groot, J. H. B., Semin, G. R. & Smeets, M. A. M. On the communicative function of body odors: A theoretical integration and review. Perspect. Psychol. Sci. 12, 306–324 (2017).Article 
PubMed 

Google Scholar 
de Groot, J. H. B. et al. A sniff of happiness. Psychol. Sci. 26, 684–700 (2015).Article 
PubMed 

Google Scholar 
Destrez, A. et al. Male mice and cows perceive human emotional chemosignals: A preliminary study. Anim. Cogn. 24, 1205–1214 (2021).Article 
PubMed 

Google Scholar 
Wilson, Id. Dogs can discriminate between human baseline and psychological stress condition odours. PLoS ONE 17, e0274143 (2022).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
D’Aniello, B., Semin, G. R., Alterisio, A., Aria, M. & Scandurra, A. Interspecies transmission of emotional information via chemosignals: from humans to dogs (Canis lupus familiaris). Anim. Cogn. 21, 67–78 (2018).Article 
PubMed 

Google Scholar 
D’Aniello, B. et al. Sex differences in the behavioral responses of dogs exposed to human chemosignals of fear and happiness. Anim. Cogn. 24, 299–309 (2021).Article 
PubMed 
PubMed Central 

Google Scholar 
Siniscalchi, M., d’Ingeo, S., Quaranta, A., D’Ingeo, S. & Quaranta, A. Lateralized emotional functioning in domestic animals. Appl. Anim. Behav. Sci. 237, 105282 (2021).Article 

Google Scholar 
Rogers, L. & Vallortigara, G. Lateralized Brain Functions: Methods in Human and Non-Human Species. vol. 122 (2017).D’Ingeo, S. et al. Horses associate individual human voices with the valence of past interactions: A behavioural and electrophysiological study. Sci. Rep. 9, 11568 (2019).Article 
ADS 
PubMed 
PubMed Central 

Google Scholar 
Siniscalchi, M., Padalino, B., Aubé, L. & Quaranta, A. Right-nostril use during sniffing at arousing stimuli produces higher cardiac activity in jumper horses. Laterality 20, (2015).De Boyer Des Roches, A., Richard-Yris, M.-A. A., Henry, S., Ezzaouïa, M. & Hausberger, M. Laterality and emotions: Visual laterality in the domestic horse (Equus caballus) differs with objects’ emotional value. Physiol. Behav. 94, 487–490 (2008).Albrecht, J. et al. Smelling chemosensory signals of males in anxious versus nonanxious condition increases state anxiety of female subjects. Chem. Senses 36, 19–27 (2011).Article 
CAS 
PubMed 

Google Scholar 
Derrickson, S. Sinister (VF)—YouTube. (Wild Bunch SA, 2001).Friard, O. & Gamba, M. BORIS: a free, versatile open-source event-logging software for video/audio coding and live observations. Methods Ecol. Evol. 7, 1325–1330 (2016).Article 

Google Scholar 
R Core Team. R: A language and environment for statistical computing (2021).Wickham, H. Ggplot2: Elegant graphics for data analysis. (2016).Hothorn, T., Winell, H., Hornik, K., van de Wiel, M. A. & Zeileis, A. coin: Conditional inference procedures in a permutation test framework (2021).Brooks, M. E. et al. glmmTMB balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling. R J. 9, 378–400 (2017).Article 

Google Scholar 
Hartig, F. DHARMa: Residual diagnostics for hierarchical (multi-level/mixed) regression models (2021).Lenth, R. V. emmeans: Estimated Marginal Means, aka Least-Squares Means. (2022).Hothersall, B., Harris, P., Sörtoft, L. & Nicol, C. J. Discrimination between conspecific odour samples in the horse (Equus caballus). Appl. Anim. Behav. Sci. 126, 37–44 (2010).Article 

Google Scholar 
Smeets, M. A. M. et al. Chemical fingerprints of emotional body odor. Metabolites 10, (2020).Sabiniewicz, A. et al. A preliminary investigation of interspecific chemosensory communication of emotions: Can Humans (Homo sapiens) recognise fear- and non-fear body odour from horses (Equus ferus caballus). Animal 11, 3499 (2021).Article 

Google Scholar 
Zhou, W. & Chen, D. Entangled chemosensory emotion and identity: Familiarity enhances detection of chemosensorily encoded emotion. Soc. Neurosci. 6, 270–276 (2011).Article 
PubMed 

Google Scholar 
Starling, M., McLean, A. & McGreevy, P. The contribution of equitation science to minimising horse-related risks to humans. Animal 6, 15 (2016).Article 

Google Scholar 
Basile, M. et al. Socially dependent auditory laterality in domestic horses (Equus caballus). Anim. Cogn. 12, 611–619 (2009).Article 
PubMed 

Google Scholar 
Hatfield, E., Cacioppo, J. T. & Rapson, R. L. Emotional contagion. Curr. Dir. Psychol. Sci. 2, 96–100 (1993).Article 

Google Scholar 
Austin, N. P. & Rogers, L. J. Limb preferences and lateralization of aggression, reactivity and vigilance in feral horses Equus caballus. Anim. Behav. 83, 239–247 (2012).Article 

Google Scholar 
Larose, C., Richard-Yris, M.-A., Hausberger, M. & Rogers, L. J. Laterality of horses associated with emotionality in novel situations Laterality Asymmetries Body. Brain Cogn. 11, 355–367 (2006).
Google Scholar 
Farmer, K., Krüger, K., Byrne, R. W. & Marr, I. Sensory laterality in affiliative interactions in domestic horses and ponies (Equus caballus). Anim. Cogn. 21, 631–637 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Chen, D. & Haviland-Jones, J. Human olfactory communication of emotion. Percept. Mot. Skills 91, 771–781 (2000).Article 
CAS 
PubMed 

Google Scholar 
de Groot, J. H. B., Semin, G. R. & Smeets, M. A. M. Chemical communication of fear: A case of male-female asymmetry. J. Exp. Psychol. Gen. 143, 1515–1525 (2014).Article 
PubMed 

Google Scholar 
Marinier, S. L., Alexander, A. J. & Waring, G. H. Flehmen behaviour in the domestic horse: Discrimination of conspecific odours. Appl. Anim. Behav. Sci. 19, 227–237 (1988).Article 

Google Scholar 
Lansade, L., Pichard, G. & Leconte, M. Sensory sensitivities: Components of a horse’s temperament dimension. Appl. Anim. Behav. Sci. 114, 534–553 (2008).Article 

Google Scholar 
Lansade, L., Bouissou, M. F. & Erhard, H. W. Fearfulness in horses: A temperament trait stable across time and situations. Appl. Anim. Behav. Sci. 115, 182–200 (2008).Article 

Google Scholar 
Hoenen, M., Wolf, O. T. & Pause, B. M. The impact of stress on odor perception. https://doi.org/10.1177/030100661668870746,366-376 (2017).Article 
PubMed 

Google Scholar 
Rørvang, M. V., Nicova, K. & Yngvesson, J. Horse odor exploration behavior is influenced by pregnancy and age. Front. Behav. Neurosci. 16, 295 (2022).Article 

Google Scholar 
Doty, R. L. & Cameron, E. L. Sex differences and reproductive hormone influences on human odor perception. Physiol. Behav. 97, 213–228 (2009).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar  More

Martins, L. C. et al. Poluição atmosférica e atendimentos por pneumonia e gripe em São Paulo, Brasil. Revista de Saúde Pública 36, 88–94 (2002).Article 
PubMed 

Google Scholar 
Goudarzi, G. et al. Health risk assessment on human exposed to heavy metals in the ambient air PM10 in Ahvaz, Southwest Iran. Int. J. Biometeorol. 62, 1075–1083 (2018).Article 
ADS 
PubMed 

Google Scholar 
Makri, A. & Stilianakis, N. I. Vulnerability to air pollution health effects. Int. J. Hygiene Environ. Health 211, 326–336 (2008).Article 

Google Scholar 
Idani, E. et al. Characteristics, sources, and health risks of atmospheric PM10-bound heavy metals in a populated Middle Eastern City. Toxin Rev. 39, 266–274 (2020).Article 

Google Scholar 
Wang, J., Hu, Z., Chen, Y., Chen, Z. & Xu, S. Contamination characteristics and possible sources of PM10 and PM2.5 in different functional areas of Shanghai, China. Atmos. Environ. 68, 221–229 (2013).Article 
ADS 
CAS 

Google Scholar 
Guarnieri, M. & Balmes, J. R. Outdoor air pollution and asthma. Lancet 383, 1581–1592 (2014).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Anderson, J. O., Thundiyil, J. G. & Stolbach, A. Clearing the air: A review of the effects of particulate matter air pollution on human health. J. Med. Toxicol. 8, 166–175 (2012).Article 
CAS 
PubMed 

Google Scholar 
Roy, D., Seo, Y.-C., Kim, S. & Oh, J. Human health risks assessment for airborne PM10-bound metals in Seoul, Korea. Environ. Sci. Pollut. Res. 26, 24247–24261 (2019).Article 
CAS 

Google Scholar 
Maesano, C. et al. Impacts on human mortality due to reductions in PM10 concentrations through different traffic scenarios in Paris, France. Sci. The Total. Environ. 698, 134257 (2020).Article 
ADS 
CAS 

Google Scholar 
Maleki, H., Sorooshian, A., Goudarzi, G., Nikfal, A. & Baneshi, M. M. Temporal profile of PM10 and associated health effects in one of the most polluted cities of the world (Ahvaz, Iran) between 2009 and 2014. Aeolian Res. 22, 135–140 (2016).Article 
ADS 
PubMed 
PubMed Central 

Google Scholar 
Medina, S., Le Tertre, A. & Saklad, M. The Apheis project: Air pollution and health—A European information system. Air Qual. Atmos. Heal. 2, 185–198 (2009).Article 

Google Scholar 
Medina, S., Plasencia, A., Ballester, F., Mücke, H. & Schwartz, J. Apheis: Public health impact of PM10 in 19 European cities. J. Epidemiol. Community Heal. 58, 831–836 (2004).Article 
CAS 

Google Scholar 
Pérez-Martínez, P. J., de Fátima Andrade, M. & de Miranda, R. M. Traffic-related air quality trends in São Paulo, Brazil. J. Geophys. Res. Atmos. 120, 6290–6304 (2015).Article 
ADS 

Google Scholar 
Sánchez-Ccoyllo, O. R. et al. Vehicular particulate matter emissions in road tunnels in Sao Paulo, Brazil. Environ. Monitoring Assess. 149, 241–249 (2009).Article 

Google Scholar 
Ribeiro, H. & de Assunção, J. V. Historical overview of air pollution in São Paulo Metropolitan Area, Brazil: Influence of mobile sources and related health effects. WIT Trans. Built Environ. 52,10 (2001).Bravo, M. A. & Bell, M. L. Spatial heterogeneity of PM10 and O3 in São Paulo, Brazil, and implications for human health studies. J. Air Waste Manag. Assoc. 61, 69–77 (2011).Article 
CAS 
PubMed 

Google Scholar 
De Freitas, E. D., Martins, L. D., da Silva Dias, P. L. & de Fátima Andrade, M. A simple photochemical module implemented in rams for tropospheric ozone concentration forecast in the metropolitan area of Sao Paulo, Brazil: Coupling and validation. Atmos. Environ. 39, 6352–6361 (2005).Article 
ADS 

Google Scholar 
Encalada-Malca, A. A., Cochachi-Bustamante, J. D., Rodrigues, P. C., Salas, R. & López-Gonzales, J. L. A spatio-temporal visualization approach of PM10 concentration data in Metropolitan Lima. Atmosphere 12, 609 (2021).Article 
ADS 

Google Scholar 
do Meio Ambiente, C. N. Institutes the national air quality control programee. Tech. Rep., Official Journal of the Federative Republic of Brazil (1989).do Meio Ambiente, C. N. Sets standards of primary and secondary air quality and even the criteria for acute episodes of air pollution. Tech. Rep., Official Journal of the Federative Republic of Brazil (1990).Artaxo, P. O estado da qualidade do ar no brasil. Work. Pap. WRI Brasil 32 (2021).Costa, A. F., Hoek, G., Brunekreef, B. & Ponce de Leon, A. C. Air pollution and deaths among elderly residents of Sao Paulo, Brazil: An analysis of mortality displacement. Environ. Health Perspectives 125, 349–354 (2017).Article 
CAS 

Google Scholar 
Bravo, M. A., Son, J., De Freitas, C. U., Gouveia, N. & Bell, M. L. Air pollution and mortality in São Paulo, Brazil: Effects of multiple pollutants and analysis of susceptible populations. J. Exposure Sci. Environ. Epidemiol. 26, 150–161 (2016).Article 
CAS 

Google Scholar 
Chiarelli, P. S. et al. The association between air pollution and blood pressure in traffic controllers in Santo André, São Paulo, Brazil. Environ. Res. 111, 650–655 (2011).Article 
CAS 

Google Scholar 
Ventura, L. M. B., de Oliveira Pinto, F., Soares, L. M., Luna, A. S. & Gioda, A. Forecast of daily PM2.5 concentrations applying artificial neural networks and holt-winters models. Air Qual. Atmos. Heal. 12, 317–325 (2019).Article 
CAS 

Google Scholar 
Leão, M. L. P., Zhang, L. & da Silva Júnior, F. M. R. Effect of particulate matter (PM2.5 and PM10) on health indicators: Climate change scenarios in a Brazilian Metropolis. Environ. Geochem. Heal. 44, 1–12 (2022).Habermann, M. & Gouveia, N. Application of land use regression to predict the concentration of inhalable particular matter in São Paulo City, Brazil. Engenharia Sanit. e Ambiental 17, 155–162 (2012).Article 

Google Scholar 
Braga, A. L. F., Pereira, L. A. A., Procópio, M., André, P. A. D. & Saldiva, P. H. D. N. Association between air pollution and respiratory and cardiovascular diseases in Itabira, Minas Gerais State. Brazil. Cadernos de Saúde Pública 23, S570–S578 (2007).Article 
PubMed 

Google Scholar 
Pinto, W. D. P., Reisen, V. A. & Monte, E. Z. Previsão da concentração de material particulado inalável, na região da grande vitória, ES, Brasil, utilizando o modelo sarimax. Engenharia Sanitária e Ambiental 23, 307–318 (2018).Article 

Google Scholar 
Schornobay-Lui, E. et al. Prediction of short and medium term PM10 concentration using artificial neural networks. Manag. Environ. Qual. An Int. J. 30, 414–436 (2018).Neto, P. S. D. M. et al. Neural-based ensembles for particulate matter forecasting. IEEE Access 9, 14470–14490 (2021).Article 

Google Scholar 
Albuquerque Filho, F. S. D., Madeiro, F., Fernandes, S. M., de Mattos Neto, P. S. & Ferreira, T. A. Time-series forecasting of pollutant concentration levels using particle swarm optimization and artificial neural networks. Química Nova 36, 783–789 (2013).Lei, T. M., Siu, S. W., Monjardino, J., Mendes, L. & Ferreira, F. Using machine learning methods to forecast air quality: A case study in Macao. Atmosphere 13, 1412 (2022).Article 
ADS 
CAS 

Google Scholar 
Yu, T. et al. Study on the regional prediction model of PM2.5 concentrations based on multi-source observations. Atmos. Pollut. Res. 13, 101363 (2022).Article 
CAS 

Google Scholar 
Li, J., Xu, G. & Cheng, X. Combining spatial pyramid pooling and long short-term memory network to predict PM2.5 concentration. Atmos. Pollut. Res. 13, 101309 (2022).Article 
CAS 

Google Scholar 
Cordova, C. H. et al. Air quality assessment and pollution forecasting using artificial neural networks in Metropolitan Lima-Peru. Sci. Rep. 11, 1–19 (2021).Article 

Google Scholar 
Plocoste, T., Calif, R. & Jacoby-Koaly, S. Temporal multiscaling characteristics of particulate matter PM10 and ground-level ozone O3 concentrations in caribbean region. Atmos. Environ. 169, 22–35 (2017).Article 
ADS 
CAS 

Google Scholar 
Calif, R. & Schmitt, F. G. Multiscaling and joint multiscaling description of the atmospheric wind speed and the aggregate power output from a wind farm. Nonlinear Process. Geophys. 21, 379–392 (2014).Article 
ADS 

Google Scholar 
Hyndman, R. J. & Khandakar, Y. Automatic time series forecasting: The forecast package for r. J. Stat. Softw. 27, 1–22 (2008).Article 

Google Scholar 
Harvey, A. C. Forecasting, structural time series models and the Kalman filter (Cambridge University Press, 1990).Book 
MATH 

Google Scholar 
Zhang, G. P. Time series forecasting using a hybrid arima and neural network model. Neurocomputing 50, 159–175 (2003).Article 
MATH 

Google Scholar 
Liao, T. W. Clustering of time series data—A survey. Pattern Recognit. 38, 1857–1874 (2005).Article 
ADS 
MATH 

Google Scholar 
Bell, M. L., Samet, J. M. & Dominici, F. Time-series studies of particulate matter. Annu. Rev. Public Heal. 25, 247–280 (2004).Article 

Google Scholar 
Hyndman, R. J. & Athanasopoulos, G. Forecasting: Principles and Practice (OTexts, 2018).
Google Scholar 
Box, G. E., Hillmer, S. C. & Tiao, G. C. Analysis and modeling of seasonal time series. in Seasonal analysis of economic time series, 309–344 (NBER, 1978).Sulandari, W., Suhartono, Subanar & Rodrigues, P. C. Exponential smoothing on modeling and forecasting multiple seasonal time series: An overview. Fluctuation Noise Lett. 20, 2130003 (2021).Article 
ADS 

Google Scholar 
Rodrigues, P. C., Awe, O. O., Pimentel, J. S. & Mahmoudvand, R. Modelling the behaviour of currency exchange rates with singular spectrum analysis and artificial neural networks. Stats 3, 137–157 (2020).Article 

Google Scholar 
Sako, K., Mpinda, B. N. & Rodrigues, P. C. Neural networks for financial time series forecasting. Entropy 24, 657 (2022).Article 
ADS 
PubMed 
PubMed Central 

Google Scholar 
Coelho, Leite et al. Statistical and artificial neural networks models for electricity consumption forecasting in the Brazilian industrial sector. Energies 15, 588 (2022).Article 

Google Scholar 
Sulandari, W., Subanar, S., Lee, M. H. & Rodrigues, P. C. Time series forecasting using singular spectrum analysis, fuzzy systems and neural networks. MethodsX 7, 101015 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
Sulandari, W. et al. Indonesian electricity load forecasting using singular spectrum analysis, fuzzy systems and neural networks. Energy 190, 116408 (2020).Article 

Google Scholar 
Rodrigues, P. C. & Mahmoudvand, R. The benefits of multivariate singular spectrum analysis over the univariate version. J. Frankl. Inst. 355, 544–564 (2018).Article 
MathSciNet 
MATH 

Google Scholar  More

Study populations and sequencing strategyDNA libraries were prepared for 1261 D. sylvestris individuals from 115 populations (5–20 individuals per population) under a modified protocol49 of the Illumina Nextera DNA library preparation kit (Supplementary Methods S1.1, Supplementary Data 1). Individuals were indexed with unique dual-indexes (IDT Illumina Nextera 10nt UDI – 384 set) from Integrated DNA Technologies Co, to avoid index-hopping50. Libraries were sequenced (150 bp paired-end sequencing) in four lanes of an Illumina NovaSeq 6000 machine at Novogene Co. This resulted in an average coverage of ca. 2x per individual. Sequenced individuals were trimmed for adapter sequences (Trimmomatic version 0.3551), mapped (BWA-MEM version 0.7.1752,53) against a reference assembly54 (ca. 440 Mb), had duplicates marked and removed (Picard Toolkit version 2.0.1; http://broadinstitute.github.io/picard), locally realigned around indels (GATK version 3.555), recalibrated for base quality scores (ATLAS version 0.956) and had overlapping read pairs clipped (bamUtil version 1.0.1457) (Supplementary Methods S1.1). Population genetic analyses were performed on the resultant BAM files via genotype likelihoods (ANGSD version 0.93358 and ATLAS versions 0.9–1.056), to accommodate the propagation of uncertainty from the raw sequence data to population genetic inference.Population genetic structure and biogeographic barriersTo investigate the genetic structure of our samples (Fig. 2A, Supplementary Fig. S2), we performed principal component analyses (PCA) on all 1261 samples (“full” dataset) via PCAngsd version 0.9859, following conversion of the mapped sequence data to ANGSD genotype likelihoods in Beagle format (Supplementary Methods S1.2). To visualise PCA results in space (Supplementary Fig. S4), individuals’ principal components were projected on a map, spatially interpolated (linear interpolation, akima R package version 0.6.260) and had the first two principal components represented as green and blue colour channels. Given that uneven sampling can bias the inference of structure in PCA, PCA was also performed on a balanced dataset comprising a common, down-sampled size of 125 individuals per geographic region (“balanced” dataset; Fig. 2B, Supplementary Fig. S3; Supplementary Methods S1.2; Supplementary Data 1). Individual admixture proportions and ancestral allele frequencies were estimated using PCAngsd (-admix model) for K = 2–6, using the balanced dataset to avoid potential biases related to imbalanced sampling22,23 and an automatic search for the optimal sparseness regularisation parameter (alpha) soft-capped to 10,000 (Supplementary Methods S1.2). To visualise ancestry proportions in space, population ancestry proportions were spatially interpolated (kriging) via code modified from Ref. 61 (Supplementary Fig. S5).To test if between-lineage admixture underlies admixture patterns inferred by PCAngsd or if the data is better explained by alternative scenarios such as recent bottlenecks, we used chromosome painting and patterns of allele sharing to construct painting palettes via the programmes MixPainter and badMIXTURE (unlinked model)28 and compared this to the PCAngsd-inferred palettes (Fig. 2B, C; Supplementary Methods S1.2). We referred to patterns of residuals between these palettes to inform of the most likely underlying demographic scenario. For assessing Alpine–Balkan palette residuals (and hence admixture), 65 individuals each from the French Alps (inferred as pure Alpine ancestry in PCAngsd), Monte Baldo (inferred with both Alpine and Balkan ancestries in PCAngsd) and Julian Alps (inferred as pure Balkan ancestry in PCAngsd) were analysed under K = 2 in PCAngsd and badMIXTURE (Fig. 2C). For assessing Apennine–Balkan admixture, 22 individuals each from the French pre-Alps (inferred as pure Apennine ancestry in PCAngsd), Tuscany (inferred with both Apennine and Balkan ancestries in PCAngsd) and Julian Alps (inferred as pure Balkan ancestry in PCAngsd) were analysed under K = 2 in PCAngsd and badMIXTURE.To construct a genetic distance tree (Supplementary Fig. S1), we first calculated pairwise genetic distances between 549 individuals (5 individuals per population for all populations) using ATLAS, employing a distance measure (weight) reflective of the number of alleles differing between the genotypes (Supplementary Methods S1.2; Supplementary Data 1). A tree was constructed from the resultant distance matrix via an initial topology defined by the BioNJ algorithm with subsequent topological moves performed via Subtree Pruning and Regrafting (SPR) in FastME version 2.1.6.162. This matrix of pairwise genetic distances was also used as input for analyses of effective migration and effective diversity surfaces in EEMS25. EEMS was run setting the number of modelled demes to 1000 (Fig. 2A, Supplementary Fig. S8). For each case, ten independent Markov chain Monte Carlo (MCMC) chains comprising 5 million iterations each were run, with a 1 million iteration burn-in, retaining every 10,000th iteration. Biogeographic barriers (Fig. 2A, Supplementary Fig. S7) were further identified via applying Monmonier’s algorithm24 on a valuated graph constructed via Delauney triangulation of population geographic coordinates, with edge values reflecting population pairwise FST; via the adegenet R package version 2.1.163. FST between all population pairs were calculated via ANGSD, employing a common sample size of 5 individuals per population (Supplementary Fig. S6; Supplementary Methods S1.2; Supplementary Data 1). 100 bootstrap runs were performed to generate a heatmap of genetic boundaries in space, from which a weighted mean line was drawn (Supplementary Fig. S7). All analyses in ANGSD were performed with the GATK (-GL 2) model, as we noticed irregularities in the site frequency spectra (SFS) with the SAMtools (-GL 1) model similar to that reported in Ref. 58 with particular BAM files. All analyses described above were performed on the full genome.Ancestral sequence reconstructionTo acquire ancestral states and polarise site-frequency spectra for use in the directionality index ψ and demographic inference, we reconstructed ancestral genome sequences at each node of the phylogenetic tree of 9 Dianthus species: D. carthusianorum, D. deltoides, D. glacialis, D. sylvestris (Apennine lineage), D. lusitanus, D. pungens, D. superbus alpestris, D. superbus superbus, and D. sylvestris (Alpine lineage). This tree topology was extracted from a detailed reconstruction of Dianthus phylogeny based on 30 taxa by Fior et al. (Fior, Luqman, Scharmann, Zemp, Zoller, Pålsson, Gargano, Wegmann & Widmer; paper in preparation) (Supplementary Methods S1.3). For ancestral sequence reconstruction, one individual per species was sequenced at medium coverage (ca. 10x), trimmed (Trimmomatic), mapped against the D. sylvestris reference assembly (BWA-MEM) and had overlapping read pairs clipped (bamUtil) (Supplementary Methods S1.3). For each species, we then generated a species-specific FASTA using GATK FastaAlternateReferenceMaker. This was achieved by replacing the reference bases at polymorphic sites with species-specific variants as identified by freebayes64 (version 1.3.1; default parameters), while masking (i.e., setting as “N”) sites (i) with zero depth and (ii) that didn’t pass the applied variant filtering criteria (i.e., that are not confidently called as polymorphic; Supplementary Methods S1.3). Species FASTA files were then combined into a multi-sample FASTA. Using this, we probabilistically reconstructed ancestral sequences at each node of the tree via PHAST (version 1.4) prequel65, using a tree model produced by PHAST phylofit under a REV substitution model and the specified tree topology (Supplementary Methods S1.3). Ancestral sequence FASTA files were then generated from the prequel results using a custom script.Expansion signalTo calculate the population pairwise directionality index ψ for the Alpine lineage, we utilised equation 1b from Peter and Slatkin (2013)31, which defines ψ in terms of the two-population site frequency spectrum (2D-SFS) (Supplementary Methods S1.4). 2D-SFS between all population pairs (10 individuals per population; Supplementary Data 1) were estimated via ANGSD and realSFS66 (Supplementary Methods S1.4), for unfolded spectra. Unfolding of spectra was achieved via polarisation with respect to the ancestral state of sites defined at the D. sylvestris (Apennine lineage) – D. sylvestris (Alpine lineage) ancestral node. Correlation of pairwise ψ and (great-circle) distance matrices was tested via a Mantel test (10,000 permutations). To infer the geographic origin of the expansion (Fig. 3), we employed a time difference of arrival (TDOA) algorithm following Peter and Slatkin (2013);31 performed via the rangeExpansion R package version 0.0.0.900031,67. We further estimated the strength of the founder of this expansion using the same package.Demographic inferenceTo evaluate the demographic history of D. sylvestris, a set of candidate demographic models was formulated. To constrain the topology of tested models, we first inferred the phylogenetic tree of the three identified evolutionary lineages of D. sylvestris (Alpine, Apennine and Balkan) as embedded within the larger phylogeny of the Eurasian Dianthus clade (note that the phylogeny from Fior et al. (Fior, Luqman, Scharmann, Zemp, Zoller, Pålsson, Gargano, Wegmann & Widmer; paper in preparation) excludes Balkan representatives of D. sylvestris). Trees were inferred based on low-coverage whole-genome sequence data of 1–2 representatives from each D. sylvestris lineage, together with whole-genome sequence data of 7 other Dianthus species, namely D. carthusianorum, D. deltoides, D. glacialis, D. lusitanus, D. pungens, D. superbus alpestris and D. superbus superbus, that were used to root the D. sylvestris clade (Supplementary Methods S1.5). We estimated distance-based phylogenies using ngsDist68 that accommodates genotype likelihoods in the estimation of genetic distances (Supplementary Methods S1.5). Genetic distances were calculated via two approaches: (i) genome-wide and (ii) along 10 kb windows. For the former, 110 bootstrap replicates were calculated by re-sampling over similar-sized genomic blocks. For the alternative strategy based on 10 kb windows, window trees were combined using ASTRAL-III version 5.6.369 to generate a genome-wide consensus tree accounting for potential gene tree discordance (Supplementary Methods S1.5). Trees were constructed from matrices of genetic distances from initial topologies defined by the BioNJ algorithm with subsequent topological moves performed via Subtree Pruning and Regrafting (SPR) in FastME version 2.1.6.162. We rooted all resultant phylogenetic trees with D. deltoides as the outgroup70. Both approaches recovered a topology with the Balkan lineage diverging prior to the Apennine and Alpine lineages (Supplementary Fig. S9). This taxon topology for D. sylvestris was supported by high ASTRAL-III posterior probabilities ( >99%), ASTRAL-III quartet scores ( >0.5) and bootstrap values ( >99%). Topologies deeper in the tree were less well-resolved (with quartet scores More