ARTICLE https://doi.org/10.1038/s43247-023-01175-9 OPEN High temperature sensitivity of monoterpene emissions from global vegetation Efstratios Bourtsoukidis 1 ✉, Andrea Pozzer 1,2, Jonathan Williams 1,2, David Makowski 3, Josep Peñuelas 4,5, Vasileios N. Matthaios6,7, Georgia Lazoglou1, Ana Maria Yañez-Serrano8, Jos Lelieveld 1,2, Philippe Ciais1,9, Mihalis Vrekoussis 1,10,11, Nikos Daskalakis 10 & Jean Sciare1 Terrestrial vegetation emits vast amounts of monoterpenes into the atmosphere, influencing ecological interactions and atmospheric chemistry. Global emissions are simulated as a function of temperature with a fixed exponential relationship (β coefficient) across forest 1234567890():,; ecosystems and environmental conditions. We applied meta-analysis algorithms on 40 years of published monoterpene emission data and show that relationship between emissions and temperature is more sensitive and intricate than previously thought. Considering the entire dataset, a higher temperature sensitivity (β = 0.13 ± 0.01 °C−1) is derived but with a linear increase with the reported coefficients of determination (R2), indicating that co-occurring environmental factors modify the temperature sensitivity of the emissions that is primarily related to the specific plant functional type (PFT). Implementing a PFT-dependent β in a biogenic emission model, coupled with a chemistry – climate model, demonstrated that atmospheric processes are exceptionally dependent on monoterpene emissions which are subject to amplified variations under rising temperatures. 1 Climate and Atmosphere Research Center (CARE-C), The Cyprus Institute, Nicosia, Cyprus. 2 Max Planck Institute for Chemistry (MPIC), Mainz, Germany. 3 French National Institute for Agriculture, Food, and Environment (INRAE), Unversity Paris-Saclay, AgroParisTech (UMR MIA), Paris, France. 4 Center for Research Ecology and Forestry Applications (CREAF), Bellaterra, Barcelona, Spain. 5 Global Ecology Unit CREAF-CSIC-UAB, Bellaterra, Barcelona, Spain. 6 Department of Public Health Policy and Systems, University of Liverpool, L69 3GB Liverpool, UK. 7 Department of Environmental Health, Harvard T.H. Chan School of Public Health, Boston, MA, USA. 8 Institute of Environmental Assessment and Water Research (IDAEA), Barcelona, Spain. 9 Laboratoire des Sciences du Climat et de l’Environnement (LSCE), Gif-sur-Yvette, France. 10 Institute of Environmental Physics (IUP), University of Bremen, Bremen, Germany. 11 Center of Marine Environmental Sciences (MARUM), University of Bremen, Bremen, Germany. ✉email: e.bourtsoukidis@cyi.ac.cy COMMUNICATIONS EARTH & ENVIRONMENT | (2024)5:23 | https://doi.org/10.1038/s43247-023-01175-9 | www.nature.com/commsenv 1 ARTICLE COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-023-01175-9 G lobal vegetation is considerably affected by rising tem- Results peratures, yet our understanding of plant response to Meta-analysis of monoterpene observations. Screening of all warming and its implications for the biosphere- indexed research articles and data sets published during 1980- atmosphere interactions remains incomplete. Biogenic volatile 2020 yielded a collection of 696 β coefficients from regression fits organic compounds (bVOCs), particularly monoterpenes (MT; obtained under diverse locations (Fig.1). The reported values C10H16), are known to be emitted in response to abiotic drivers varied widely, ranging from −0.07 to 2.39 °C−1. To identify the such as temperature1–3. Due to their essential role in warming- factors that explain this variability, we collected 35 additional induced responses and emission composition characteristics, MTs parameters that described the experimental procedure, location, have even been attributed the name ‘thermometer of plants’4. time, vegetation state (i.e. the plant functional type and tree Interlinked to their temperature-dependent emission, the multi- characteristics) and the characteristics of MT emission, including faceted roles of MTs in the earth’s system extend from commu- the conclusion of each study (Supplementary Table 1). nication signals to influencing atmospheric processes5,6. Almost all observations were in the Northern Hemisphere, and Once released into the atmosphere, they become key players in the majority of the studies were conducted under field conditions, atmospheric chemistry and physics7,8. Due to their high reactivity using branch enclosures and offline chemical analyses for with atmospheric radicals, they are key contributors to the for- quantifying VOCs using adsorbent tubes (Supplementary Fig. 1). mation and growth of secondary organic aerosols (SOA), indir- The derivation of the rate of emission (enclosures vs eddy ectly impacting the radiative balance of the atmosphere9–11. In covariance) and the techniques of sampling and analyzing VOCs the presence of nitrogen oxides, MTs also participate in the for- did not differ β amongst the studies (offline vs. online), validating mation of tropospheric ozone, a harmful atmospheric pollutant12. the upscaling of β from leaf to canopy levels (Supplementary Their role on atmospheric oxidation and SOA formation make Fig. 2). In the majority of cases, the reported MT emissions were their accurate simulation essential for understanding and pre- solely temperature-dependent, accounting for 499 instances or dicting the impacts of climate change. 72% of the total reports. Emissions dependent on both Global emission models employ empirical algorithms to temperature and light were reported in 142 instances (20.4%), simulate variations in MT emission rates13. These algorithms while the remaining 55 β-coefficients (7.9%) were derived when encompass a temperature response mechanism rooted in enzy- authors identified both temperature and another environmental matic activity and a light response mechanism based on electron driver as influencing the emissions (Supplementary Fig. 3). The transport14. Together, they offer a straightforward yet mechan- three distinct clusters exhibited the following median values: istic approach to MT emission modeling. Generally, rising tem- temperature alone (0.11 ± 0.13 °C−1), temperature and light peratures stimulate the emission of MTs stored in leaf pools, (0.14 ± 0.12 °C−1), and temperature combined with another while both temperature and light govern the de novo production environmental driver (0.10 ± 0.05 °C−1) (Supplementary Fig. 4a). and subsequent release of these compounds15. Recent scientific The additional environmental drivers reported include relative progress has prompted the incorporation of additional para- humidity, soil water content, CO2 levels, seasonality, transpira- meters into the model, such as leaf age, soil moisture, leaf area tion, and photosynthesis. The clearest observed pattern was the index, and CO2 inhibition (refer to Supplementary information). increase in β with the coefficient of determination (R2) obtained Despite these advancements, the majority of studies still attributes from the regression fits between emission and temperature for temperature as the primary driver of MT emissions from global each experimental dataset. In this context, R2 serves as a statistical vegetation, commonly represented by the following equation:13: measure that characterizes the goodness of the fit of Eq. (1), elucidating the degree to which the parameterization explains the  intrinsic association between temperature and emission rates. EMT ¼ E30 exp β T  T30 ð1Þ Our meta-analysis, based on a standard random-effect statistical model19, determined the values of β within 95% confidence intervals, across R2 bins (Fig. 2). A linear relationship was Here, E30 is the emission potential under standardized condi- identified between the quality of the regression fit and β for MT tions (30 °C), and β is an empirical coefficient (in °C−1) derived emissions that had larger values as the quality of the fit increased. from the fit of the regression between measured emissions and The majority of the data (ca. 52%) had an R2 > 0.6, and the temperature (T). It should be noted that MT emissions are best reported dependencies on temperature were significantly higher described by leaf temperature since leaf and air temperatures can than the established value of 0.10 °C−1. This observation indicates vary considerably among plant species due to the different plant- a clear underestimation of a globally uniform β. specific ability to cool down plant tissues16. E30 is a Plant Func- The linear increase of β with the goodness of fit of Eq. (1) (i.e. tional Type (PFT)-specific parameter that defines the MT emis- R2) reveals that when temperature is the dominant driver, the sions at 30 degrees Celsius (or 303.15 K), independent of the temperature sensitivity reaches its maximum. In fact, when temperature responses from the individual plants. These are considering instances with R2 > 0.9, the β-coefficient is specified by the slope of the exponent (i.e. the β-coefficient) and 0.17 ± 0.3 °C−1, 70% higher than the value currently consid- are notable indicators of the plant responses to temperature. ered. Conversely, at mid- and lower R2 values, β declines, To date, global emission models, such as MEGAN13, use a demonstrating a decrease in temperature sensitivity as the fixed dependence on temperature (β = 0.10 °C−1) for all types of emission variance becomes less explainable by temperature. vegetation and ecosystems. The purpose of this meta-analysis is to This implies that other environmental drivers regulate the conduct a comprehensive evaluation of this parameter by col- emission rates. Considering the emissions that were reported to lecting and analyzing all experimentally derived β coefficients be dependent on both temperature and light, β and R2 followed published in the four decades since the first observation of this the same trend. However, the values for the temperature and relationship in 198017. Moreover, we aim to assess potential light (T + L) cluster were consistently higher compared to the refinements of β and investigate the resulting impact on atmo- temperature-only (T) cluster (Supplementary Fig 4b). An spheric chemistry and physics by performing coupled simulations ANOVA test conducted on these two datasets yielded a between the most established emission model for biogenic trace p-value of 0.017, indicating statistically significant differences gases (MEGAN)13 and an atmospheric chemistry–climate model in temperature responses when MTs are also emitted as a (EMAC)18. function of light. 2 COMMUNICATIONS EARTH & ENVIRONMENT | (2024)5:23 | https://doi.org/10.1038/s43247-023-01175-9 | www.nature.com/commsenv COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-023-01175-9 ARTICLE Fig. 1 Locations of experimentally derived β. The size of the circles indicates the number of coefficients obtained from the literature. The studies considered can be found as a separate table in the Supplementary information. Ecosystem Functional Type Latitude Longitude Tree age E30 Elevation T (average) Year Month Sampling height Vegetation Season Location Tree height Time 0 10 20 30 40 50 60 70 80 90 100 Normalized importance in -factor Fig. 2 Groups of coefficients of determination (R2) for experimentally derived dependencies (β) of MT emissions on temperature. The blue and Fig. 3 Feature importance ranking. Indicative contribution of individual green error bars indicate 95% confidence intervals. On the right, the parameters of vegetation, location, and time parameters to the variation of distribution of the N = 696 values of β extracted from the literature is the β coefficient. The parameters are explained in Supplementary Table 1. displayed together with a boxplot that illustrates their median (white circle) and 25th and 75th percentiles (lower and upper bounds of the box), respectively. Data below the red dashed line (on the right) were obtained plant functional type (PFT) as defined by Guenther et al.13, is a from regressions with poor goodness of fit (R2 < 0.2) and were excluded critical factor that affects the temperature sensitivity of MT from further analysis. emissions. The geographic location (latitude and longitude) remained associated with the vegetation types that grow within Interestingly, β obtained from field data was systematically the corresponding geographically distributed ecosystems. These higher than β obtained under controlled conditions, suggesting results, however, suggest that the temperature responses of MT amplified sensitivities to temperature under natural conditions. emissions may exhibit dynamic variations throughout the tree’s This discrepancy was further enhanced at high R2 values, where age. Finally, the FIR analysis did not clearly identify seasonality or the MT emissions were primarily linked to temperature. Based on a relationship between tree or sampling height and the coefficient this observation, it can be concluded that laboratory studies may of temperature dependency. However, a relatively consistent ratio underestimate the temperature sensitivity in real-life conditions. between sampling and tree height has been identified (0.7 ± 0.2), To better understand the key factors influencing β values, revealing the preference towards sun-exposed parts of the tree we employed machine learning techniques, specifically the (Supplementary Fig. 5). Summarizing, it should be noted that Feature Importance Ranking (FIR) tool. This approach evaluates while a clear the relationship between PFT and the β coefficient the significance of individual parameters, such as vegetation type, became evident from the FIR analysis, the role of the other location, and time variables, in the determination of β (Fig. 3). parameters might have underappreciated due to the frequently The results suggest that the type of ecosystem, specifically the unreported values. COMMUNICATIONS EARTH & ENVIRONMENT | (2024)5:23 | https://doi.org/10.1038/s43247-023-01175-9 | www.nature.com/commsenv 3 ARTICLE COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-023-01175-9 Fig. 4 Plant functional type-dependent β coefficients. Dependence on temperature (β) for the plant functional types (a) and wider categories of woody plants (b). The error bars in panel a indicate 95% confidence intervals, and the size of the blue bullet points is proportional to the global surface area of each plant functional type. NETe needleleaf evergreen temperate forest, NEB needleleaf evergreen boreal forest, NDB needleleaf deciduous boreal forest, BETr broadleaf evergreen tropical forest, BETe broadleaf evergreen temperate forest/shrubs, BDTe broadleaf deciduous temperate forest/shrubs, BDB broadleaf deciduous boreal forest/shrubs. Revised temperature-dependent coefficients (β). The statistical BASE simulation, where MEGAN was used in its standard con- model used in this study analyzed the entire dataset consisting of figuration (β = 0.10 °C−1), 2) an ALL simulation, where MEGAN 559 data points with an R2 value greater than 0.2 (see methods). was used in its standard configuration but with β set at 0.13 °C−1, The results showed a uniform value of 0.13 ± 0.1 °C−1 for β, with 3) a MTRP simulation, where β was based on the monoterpene different temperature-dependent coefficients for each PFT, as species (Supplementary Fig. 7), and 4) a PFT simulation, where β seen in Fig. 4. Removing the lowest 10% of data points had little in MEGAN varied based on PFT (Fig. 4, Supplementary Code 1). effect on the revised values, as demonstrated in Supplementary The new simulations resulted in significant and comparable Fig. 6. The analysis showed that, except for boreal needle leaf changes in emission rates and atmospheric feedbacks compared deciduous trees and broadleaf deciduous shrubs, almost all other to the current model (Supplementary Fig. 8). Based on the results PFTs have significantly larger values than 0.10 °C−1. This is presented in Figs. 3 and 4, as well as Supplementary Figs. 6, 7, we particularly pronounced for the boreal needle leaf evergreen determined that the most appropriate strategy was to revise the forests (βΝΕΒ = 0.15 °C−1; range 0.12–0.19 °C−1) and the tropical temperature sensitivity coefficient (β) across different plant broadleaf evergreen forests (βBETr = 0.20 °C−1; range functional types (PFTs). 0.14–0.3 °C−1). Generally, broadleaf trees display higher tem- In the “PFT” simulation, global annual average MT emissions perature dependencies in the release of MTs by their foliage decreased by 13%, with the most pronounced reductions compared to needle-leaf trees and broadleaf shrubs. Notably, no β occurring for needle-leaf evergreen boreal trees (44%), broadleaf values were reported for agricultural ecosystems or grasslands, deciduous boreal trees (44%), and broadleaf evergreen temperate highlighting the need for further research to determine the trees (41%) (Fig. 5, Supplementary Table 2). The largest seasonal temperature responses for these climate-sensitive ecosystems. increase was observed in Siberia during winter, which is A β coefficient was derived for each of the reported MT species, attributed to the lower β values from broadleaf deciduous shrubs considering the importance of processes of atmospheric oxidation and generally low MT emissions in this region. It should be noted due to their different rates of reaction with atmospheric radicals. that the β applied for this PFT is grounded in findings from a Supplementary Fig. 7 shows that, although β was higher for Δ3- single study20. Although MT emissions from this ecosystem carene and trans-β-ocimene, no significant differences were remain relatively low21, emerging research indicates that the identified across all MTs. However, it is important to note that subarctic tundra might exhibit higher sensitivity to temperature the majority of studies (74%) reported the sum of MTs, which changes in both MT and isoprene emissions22–24. Given limits the ability to draw a clear conclusion on the temperature the accelerated warming of this delicate ecosystem compared to sensitivity of speciated MTs. This limitation is further exacerbated the global average25, it is crucial to gather more data from the by the fact that Proton Transfer Reaction instruments measure pan-Arctic region to accurately characterize the temperature the cumulative total of all MTs, whereas Gas Chromatography responses of regional vegetation. methods only account for selected species in their reported sums. In contrast, tropical ecosystems, especially the Amazon rain- Summarizing the reconsiderations needed for determining β forest, had higher annual β values, resulting in local increases of for the MTs, we concluded that the most sensible approach was to up to 30% in MT emissions. The simulations also showed that assign a different value for each PFT. Alternatively, a uniform during the Amazonian dry season, MT emissions could even value of 0.13 °C−1 could be considered, and potential discrepan- double despite similar annual averages (Supplementary Fig. 9). cies for different monoterpenes may not be disregarded due to Moreover, the investigation of daily extremes for each season insufficient published data. revealed occasional increases of more than three-fold in both tropical and arctic ecosystems (Supplementary Fig. 10). In light of the higher temperature sensitivity of MT emissions, it becomes Model simulations. To comprehensively evaluate the implica- evident that the increasing frequency and strength of heat waves26 tions of revised MT emission temperature dependencies, we will markedly amplify the warming-induced MT emissions. performed coupled simulations employing the most widely used The most pronounced difference between the two simulations VOC emission model (MEGAN)13 and a state-of-the-art global was the diurnal variation in the rates of emission (Supplementary atmospheric chemistry–climate model (EMAC)18. The simula- Fig. 11, 12). Tropical forests exhibited especially marked diel tions were performed at 1.8 by 1.8-degree resolution and with cycles, which could be further enhanced over seasonal and daily hourly outputs. In total, we performed four simulations: 1) a timescales. The maximum differences in daily standard deviation 4 COMMUNICATIONS EARTH & ENVIRONMENT | (2024)5:23 | https://doi.org/10.1038/s43247-023-01175-9 | www.nature.com/commsenv COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-023-01175-9 ARTICLE a b Mean relative difference of annual emissions (in %) Mean relative difference of daily standard deviations (in %) c d Mean relative difference of daytime masked OH (in %) Mean relative difference of annual SOA (in %) Fig. 5 Model simulations. Annual mean relative differences between the BASE simulation, where MEGAN was used in its standard configuration (β = 0.10 °C−1), and a PFT-dependent β, based on the median values presented in Fig. 4. The differences are averages of hourly data over a climatic year. a Total monoterpene emissions, b The standard deviation of daily emissions, c Daytime OH, and d SOA production. Fig. 6 Monoterpene emissions over temperature gradients for the three main ecosystems. a Emissions using the standardized potential emission at 30 °C and β in MEGAN v2.1 (dashed lines) and using β derived from this study (panel b and Supplementary Table 2). The maps in the right panel represent the global coverage of boreal (blue), temperate (green), and tropical (red) forests. exceeded 300%, highlighting the substantial variability in MT emissions between the model and observations in tropical emissions due to temperature dependencies. We applied Eq. (1) ecosystems28. using data for temperature collected at the Amazon Tall Tower The relationship between temperature and MT emission Observatory (ATTO) in 2014 and 201527 to demonstrate the follows an exponential function that is standardized at 30 °C dynamics of PFT-dependent β (Supplementary Fig. 13). The (Eq. (1)). Revisions of the β values according to their PFT resulted increase in β (from 0.10 to 0.20 °C−1) led to slightly lower average in substantial changes in the simulated emissions occurred at emissions because the average temperature is lower than 30°C either low or very high temperatures (Fig. 6, Supplementary (see Fig. 6). The variations in the standard deviation of the hourly Fig. 14). We derived β weighted by area coverage for each PFT to simulated emissions, however, averaged twice as large and were illustrate the simplified dynamics of emission for the three particularly pronounced under extreme conditions such as the El dominant ecosystems (Supplementary Table 3). As shown in Niño year of 2015. Such amplified diurnal variations help to Fig. 6a and Supplementary Fig. 14, the current models of explain the discrepancies in the standard deviation of the emission have clearly overestimated MT emissions under 30 °C. COMMUNICATIONS EARTH & ENVIRONMENT | (2024)5:23 | https://doi.org/10.1038/s43247-023-01175-9 | www.nature.com/commsenv 5 ARTICLE COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-023-01175-9 In contrast, the current models have underestimated the corresponding rise in the production of SOA, thereby enhancing emissions for temperatures above 30 °C, which may help to the radiative cooling effects. These findings emphasize the urgent explain the observed discrepancies in diurnal variation of MT need for accurate simulations of the temperature effects on MT emission in warm environments28,29. emissions to improve our understanding of atmospheric chemistry The parameter E30, representing the emission potential at and physics, and ultimately, to better predict the impacts of standard conditions, plays a crucial role in determining the rate of climate change. MT emissions. We chose to perform the simulations adopting the values in MEGAN because our analysis identified large uncertainties in the potential emissions, which increased further Discussion due to the uncertainties in the conversion of units (typically from As global temperatures continue to rise and extreme heat events ng C g(dry weight)−1 h−1 to μg C m−2 h−1) and the leaf area become more frequent26, the temperature sensitivity of forests is index derived from the literature and satellite observations emerging as a critical issue for understanding the impacts of (Supplementary Fig. 15, 16). climate change on forest ecosystems and atmospheric chemistry. One important aspect of this sensitivity is the emissions of bVOCs, such as monoterpenes, which will increase with rising Atmospheric implications. The direct implications on the pro- temperatures. The stronger emissions resulting from rising tem- cesses of atmospheric oxidation were evaluated by investigating peratures can have far-reaching and uncertain consequences for the differences in concentration between model simulations. The the biosphere, potentially disrupting its delicate balance and hydroxyl radical (OH) is the most important daytime oxidant feedback mechanisms between ecology, atmospheric chemistry that reacts strongly with most organic molecules, particularly and climate6. Therefore, the accurate simulation of MT emissions MTs. The model simulations showed that using a PFT-dependent is essential for evaluating plant responses and the respective β increases OH concentrations globally (Fig. 5). The boreal-forest feedbacks and improving future projections. belt was particularly affected, with ca. 12% higher annual average To date, various algorithms have been developed to simulate OH concentrations compared to the BASE simulation MT emissions, which account for either short-term volatilization (β = 0.10 °C−1). Increased OH concentrations in NH summer (emission-based)13,17 or long-term production of MT linked to were particularly pronounced (up to 20%) (Supplementary photosynthesis (production-based)32,33. Emission-based models, Fig. 17). Substantial daytime changes in MT emissions in EMAC such as MEGAN, are commonly used to describe observations at indicated that OH concentrations could occasionally be even as different scales, from leaves and branches to canopies, by simu- much as double compared to what is currently simulated (Sup- lating global MT emissions from terrestrial vegetation. These plementary Fig. 18). It is therefore concluded that increasing the models account for multiple parameters, including temperature, temperature sensitivity of MT emissions will result in a global light, foliar age, soil moisture, and CO2 inhibition. In this meta- increase of the atmospheric oxidative capacity of OH, demon- analysis, we focused on the effect of temperature, particularly in strating the important role of MTs in atmospheric chemistry. The evaluating the warming emission-response that defines tem- implications of an enhanced OH abundance over the northern perature sensitivity (β) of MT emissions. latitudes include a reduction in the atmospheric lifetime of The β coefficient for the temperature-dependent term in methane, which is prolifically emitted in these regions30. Eq. (1) is shown to vary across seasons34,35, tree species36, tree In contrast to OH, average tropospheric concentrations of age37, developmental stage38, ambient humidity39, level of pho- ozone (O3) were only moderately affected as the maximum O3 tosynthetically active radiation40, CO2 abundance41, soil concentration change in the summer was within 4−5% for the moisture42, and other drivers43. Both the FIR analysis and sim- Northern Hemisphere (Supplementary Fig. 19). Regional changes plified correlation analyses on the new data collected in this study in O3 concentration on specific days, however, could increase by demonstrated that the vegetation type was the most important up to 10% (Supplementary Fig. 20). Monoterpenes contributed regulator of the sensitivity to temperature of MT emissions from both to the production (in the presence of NOx) and chemical vegetation. On average, plants that grow in warmer ecosystems removal of tropospheric O3 by reaction schemes that could differ appear more sensitive, indicating the adjustment of plants in considerably between locations (e.g. VOC or NOx limited response to warming and applying a PFT-dependent β helps to environments31). By reducing MT emissions in our updated explain the frequently overestimated MT emissions in several emission model, O3 loss was generally minimized, increasing its environments29,44–47. abundance globally, except for the hot and dry season of the Using a PFT-dependent β is a considerable step forward in Amazon rainforest when the high ambient MT concentrations understanding the sensitivity of MT emissions to temperature. remove O3 at higher rates compared to the base model However, given the complexity of nature, it is still necessary to simulations. acknowledge the remaining uncertainties in estimating and Changes in the dynamics of atmospheric oxidation over forests modeling the β for MT emissions. For example, our findings influenced the yield of secondary organic aerosol particles around reveal that the sensitivity of MT emissions to temperature is also the globe (Fig. 5, Supplementary Fig. 21). The revised simulations influenced by co-occurring environmental drivers and is different indicated lower SOA production over boreal and temperate for laboratory and field conditions. On the same lines, it has been forests. Seasonal averages nonetheless indicated increased SOA recently shown that the β coefficient can significantly increase formation over both tropical forests and oceans. under heat stress conditions48. Creating an accurate para- The study’s results underscore the crucial role of the meterization that considers the diverse and frequently co- temperature dependency of MT emissions in shaping atmospheric occurring drivers of emissions is a challenging task since plants oxidation processes and reaction products. The model simulations have developed distinct coping mechanisms to deal with envir- demonstrated that this sensitive coefficient plays an important role onmental stressors. It has been demonstrated that even plants of in atmospheric processes that are amplified seasonally and have a the same species and identical growing conditions can exhibit significant impact on the formation of secondary organic aerosol large chemodiversity in their MT emissions, adding to the com- particles globally. Given the higher temperature sensitivities plexity of plant responses49. presented here, the rising global temperatures will amplify MT Besides the temperature sensitivity parameter (β), other para- emissions. Consequently, this increase in MTs will lead to a meters also play a role in the emissions of MTs by terrestrial 6 COMMUNICATIONS EARTH & ENVIRONMENT | (2024)5:23 | https://doi.org/10.1038/s43247-023-01175-9 | www.nature.com/commsenv COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-023-01175-9 ARTICLE vegetation. One critical parameter is the potential emission at by extracting all available information for 35 parameters standard conditions (E30), which regulates the rate of emission by (Supplementary Table 1), vectorising them (i.e. annotated/ forests. Our meta-analysis collected E30 values from different assigned a number to each character class), and then proceeded ecosystems (see Supplementary Figs. 15, 16) and found that they with our statistical analyses. vary greatly by orders of magnitude, revealing substantial We categorized the main conclusions of each study into three uncertainties in the magnitude of simulated emissions. Given that groups based on the factors driving monoterpene emissions: both β and E30 define the emission rate of MTs from global temperature alone, temperature and light, and temperature in vegetation, further research is necessary to validate the current combination with other environmental factors such as soil approach of implementing constant E30 values in emission moisture and relative humidity. The experimental techniques models. For example, although CO2 concentrations can impact employed in these studies varied, including different environ- both β and E30 values, limited data in the reviewed literature mental conditions, methods for determining emission rates, VOC hindered a conclusive assessment. Moreover, uncertainties in the sampling, and chemical analyses (Supplementary Fig. 2). How- leaf area index and the lack of studies in widespread ecosystems, ever, these variables were comprehensively reported in all studies, including agricultural and grasslands, add to the unknowns and making the annotation process relatively straightforward. We shape the direction of future research. annotated the β data points based on the method that was used to The estimation of MT emissions has significant implications collect the majority of the data, unless there were two different for atmospheric processes, specifically for the mechanisms of sampling methods used. In such cases, we used the annotation atmospheric oxidation and subsequent SOA formation. Our method that corresponded to the primary sampling method. study demonstrates that increasing β leads to a decrease in MT Additionally, when gas chromatography was used for measuring emissions under mean temperatures below 30 °C, resulting in an individual monoterpenes, we annotated the regression fits increase in global OH concentrations. This observation is espe- accordingly. cially relevant given the existing discrepancies between measured The year, month, and season of the experiments were assessed. and simulated OH concentrations in unpolluted forested Typically, the experiments were conducted within a single year. regions50. Our updated simulations, which incorporate the PFT- In cases where the experiments were carried out over two based temperature sensitivity of MT emissions, provide an consecutive years and regression fits were applied for all available explanation for these discrepancies and underscore the strong data, we used the year in which the majority of the observations link between MT emissions and atmospheric processes. These were collected. If the experiments were conducted in the same findings not only enhance our understanding of the complex month of two different years, we annotated the year of the first mechanisms of atmospheric oxidation but also highlight the observation. For experiments where the year was not reported critical importance of accurately estimating MT emissions in (typical for laboratory experiments), we annotated the year prior simulating the behavior of atmospheric processes. to publication. The month of an experiment was annotated with With newly discovered biogenic MT sources27,51–54 and the numbers 1-12, and 13 was used for longer and/or mixed periods challenge of modeling co-occurring environmental drivers on the across years. Seasonal observations and laboratory experiments biosphere, our study highlights the need for more process- were annotated differently. The seasons were annotated accord- oriented research of biosphere-atmosphere interactions, particu- ingly in this classification. larly in tropical, pan-Arctic, grassland, and agricultural ecosys- The exact location of the field experiments included the tems. As the effects of climate change intensify, biogenic VOC latitude, longitude, and meters above sea level. If the authors emissions from global vegetation will play a crucial role in eval- reported a specific station instead of the exact coordinates, we uating the health of ecosystems and influencing the atmospheric used the cited literature to derive this information, but we used oxidation capacity, with implications for the chemical composi- Google Earth and derived the coordinates if deriving this tion, aerosols, and climate. information was not possible. Information for the tree species (both Latin and common names) was collected from each publication. Considering the geographical information, we Methods assigned each tree species to its plant functional type (PFT) as Experimental design. Searching the Web of Science, selecting defined in the Model of Emissions of Gases and Aerosols from results from the WOS, BCIBIOSIS, CCC, DRCI, and RSCI Nature (MEGAN)13. Tree age (in years), tree height (in meters), databases, and using the keywords “monoterpenes”, “emissions”, and sampling height (in meters) were rarely reported together, so and/or “temperature” identified 745 peer-reviewed studies pub- several data gaps could not be filled. Finally, the minimum, lished between 1980 and 2020. Screening of these articles (page by maximum, mean, and range of atmospheric temperatures were page) decreased the number of studies to 84 that reported β filled as reported, but if not reported, the data were collected from derived from regression fits of experimental observations. The the publication figures when possible. data available for tropical plants were insufficient for conducting These studies had different scientific objectives, so the reliable statistical analyses. To address this limitation, we included collection and classification of all 35 parameters was a strenuous three additional studies that reported relevant data. The timelines task. All parameters were cross-checked by more than one author for temperature and MT emission were extracted from published of this study to avoid mistakes during collection. plots in Jardine et. al4, and Langford et al.55 using a web-based tool (WebPlotDigitizer; https://automeris.io). Data from the third study were directly obtained from Yáñez-Serrano et al.56 We compiled a data set consisting of 696 values of β. To the best of Statistical analysis. Considering that a low R2 indicates a weak to our knowledge, we accounted for all values reported in the lit- non-existent relationship between MT emissions and tempera- erature. We may, however, have missed some, mainly because ture, we chose to disregard the lowest 10% of the reported β (56 they were not appropriately indexed in the literature data bases. values), which represented R2 < 0.2. This approach increased the The studies used in this meta-analysis are listed in the Supple- quality of our data set while retaining a substantial number of mentary Information (Supplementary References list). observations for further statistical analyses. R2 was not reported We investigated potential relationships with experimental, in 81 cases, so the overall data that we considered for analysis are geographical, plant-specific, seasonal, and regression-fit variables comprised of 559 experimental values. COMMUNICATIONS EARTH & ENVIRONMENT | (2024)5:23 | https://doi.org/10.1038/s43247-023-01175-9 | www.nature.com/commsenv 7 ARTICLE COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-023-01175-9 The remaining data were analyzed using random-effects RMSE and the greater R2 out of the three (Supplementary models, which are commonly used in meta-analyses57–59. A first Fig. 22). model without a covariate was fitted to estimate the global mean β across all studies. Another model with functional type as a categorical covariate was then fitted to the data set for analyzing Model simulations. We used the global ECHAM/MESSy the relationship between β and PFT. Finally, a categorical variable Atmospheric chemistry – Climate (EMAC) model, which simu- including categories of R2 covering intervals of 0.1 was defined lates atmospheric chemical and meteorological processes and and included as a covariate in another random-effect model for interactions with oceans and the biosphere71,72. The model uses studying the relationship between β and R2. This model was fitted the second version of the Modular Earth Submodel System to the data collected under either field or controlled conditions (MESSy2) to link multi-institutional computer codes. The core for analyzing the sensitivity of β to the type of experiment. A atmospheric model was the 5th generation European Centre/ random study effect was included in all fitted models to account Hamburg general circulation model (ECHAM5), into which for the heterogeneity between studies. The models were fitted updates and improvements in boundary layer, radiation, and using the lmer function of the lme4 R package by restricted convection routines have been introduced72–74. Additional maximum likelihood60. The accuracy of estimated β was assessed descriptions, references, and information for the model are by calculating 95% confidence intervals. All calculations were available at https://www.messy-interface.org. We applied EMAC implemented using R v4.1.2. (ECHAM5 version 5.3.02, MESSy version 2.55.0) with a spherical truncation of T63 (corresponding to a quadratic Gaussian grid of approximately 1.8 × 1.8 degrees of latitude and longitude) with 31 Unit conversions for potential emission at standardized con- levels of vertical hybrid pressure to 10 hPa. ditions (30 °C, E30). The rates of emission of monoterpenes were The various submodels represented tropospheric processes and typically reported in units of carbon mass per gramme of foliar their interactions with oceans, land, and human influences dry weight per hour (ng C g(dry weight)−1 h−1). The biogenic describing emissions, including isotopic composition, radiative emission factors in MEGAN are in μg m−2 h−1, so we used the processes, atmospheric multiphase chemistry, aerosols, and leaf area index (LAI) and the specific leaf area (SLA) for each mechanisms of deposition73,75. The set-up used in this simulation point to convert the reported values. SLA was obtained from the was identical to that used by Pozzer et al.18, where a detailed TRY database61 (https://try-db.org, last accessed 8 June 2022). scheme of the oxidation of volatile organic compounds (VOCs) The data sets used were 1-km62 maps scaled up from trait data (the Mainz Organic Mechanism) was coupled to a base set of measured in situ. The LAI data we used were the monthly cli- volatility (ORACLE76) to simulate the partitioning of organic matological data from the ORNL DAAC database63. The data for gases and aerosols in unprecedented detail for a SLA and LAI were extracted from the above data sets, for the chemistry–climate model. We only briefly summarized the most longitudes and latitudes of our data points while in case of important characteristics (see Pozzer et al.18 for more details). missing values, the nearest neighbor with valid data was used. EMAC simulates gas-phase and heterogeneous chemistry using the MECCA submodel, which accounts for the photochemical Machine learning. We used machine learning methods to oxidation of natural and anthropogenic VOCs76–78. Processes of examine the relationship and the importance of β coefficient to aerosol microphysics and gas/aerosol partitioning were simulated the three broad categories (vegetation type, location and time using the GMXe submodel79. The distribution of aerosol sizes was variables) identified in our meta-analysis. As each of these broad described using seven interacting lognormal modes (four categories had sub-divisions, we examined them using the Fea- hydrophilic and three hydrophobic modes). The composition of tured Importance Ranking (FIR) approach of machine learning. aerosols was uniform within each mode (internally mixed) but Feature (or variable) importance ranking refers to a task that could vary between modes (externally mixed). The four modes of measures the contributions of individual input features (variables) hydrophilic size (nucleation, Aitken, accumulation, and coarse) to the performance of a supervised learning model64, and effec- encompassed the spectrum of aerosol sizes. The composition of tively addresses inter-correlated parameters. Feature importance inorganic aerosols was determined using the ISORROPIA-II ranking has become one of the most powerful tools in explain- thermodynamic equilibrium submodel80, which calculates the able/interpretable models to facilitate understanding and dis- gas/liquid/solid equilibrium partitioning of inorganic compounds covery of key factors in a specific domain65,66. and water. The components of aeolian dust can exist in the form Specifically, we used Caret and Cubist packages in R with of mineral salts in the solid phase and ions in the aqueous tuning two hyper-parameters: neighbors (#Instances) and phase81. The composition and atmospheric evolution of organic committees (#Committees), which are the ones to most likely aerosol compounds were simulated using the ORACLE sub- have the largest effect on the final performance of the Cubist model, which represents classes of the volatility of organics by model67. Cubist is a rule-based tree algorithm, where a tree is their effective saturation concentrations82. The biogenic emis- grown, and the terminal leaves contain regression models. These sions of non-methane VOCs were calculated online using models are based on the predictors used in previous splits. In MEGAN13. The results of the model for the last decade have these algorithms, the prediction is made using the linear been extensively tested against measured data for gases and regression model at the terminal node of the tree but is smoothed particles from ground-based networks monitoring air quality and by taking into account the prediction from another model in the global observations from satellites83–86. previous node of the tree (which also occurs recursively up the The dynamics were not nudged, unlike the study by Pozzer tree). The tree is reduced to a set of rules, which initially are paths et al.18, but the model ran freely, forced only by climatological from the top of the tree to the bottom68. The performance is sea-surface temperature (SST) and sea-ice coverage (SIC) taken from each combination of the hyper-parameters tuning obtained from the ERA5 data for 2010–2019, which allowed the with the grid search method with cross-validation (CV)69. To model to recreate a climatological meteorology without strong avoid bias in data selection, we applied 10-fold CV69,70, while the extremes in the forcing (SST and SIC). The meteorology/radiation model’s final performance against the test dataset was validated and the chemistry were also fully decoupled, so all simulations using R2 and Root Mean Square Error (RMSE). We also tested an performed developed identical binary meteorology, allowing ANN and a random forest algorithm but Cubist had the smaller comparisons between the simulations. In summary, any 8 COMMUNICATIONS EARTH & ENVIRONMENT | (2024)5:23 | https://doi.org/10.1038/s43247-023-01175-9 | www.nature.com/commsenv COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-023-01175-9 ARTICLE difference between the concentrations of MTs in the simulations 22. Li, T., Baggesen, N., Seco, R. & Rinnan, R. Seasonal and diel patterns of were only due to changes in the emissions and not in different biogenic volatile organic compound fluxes in a subarctic tundra. Atmos. modes of transport. Environ. 292, 119430 (2023). 23. Seco, R. et al. Strong isoprene emission response to temperature in tundra vegetation. Proc. Natl. Acad. Sci. USA. 119, e2118014119 (2022). 24. Rinnan, R. et al. Separating direct and indirect effects of rising temperatures Data availability on biogenic volatile emissions in the Arctic. Proc. Natl. Acad. Sci. USA. 117, The meta-analysis dataset is openly available at https://doi.org/10.5281/zenodo. 32476–32483 (2020). 10043961. 25. Rantanen, M. et al. The Arctic has warmed nearly four times faster than the globe since 1979. Commun. Earth Environ. 3, 168 (2022). 26. Perkins-Kirkpatrick, S. E. & Lewis, S. C. Increasing trends in regional heatwaves. Nat. Commun. 11, 3357 (2020). Code availability 27. Bourtsoukidis, E. et al. Strong sesquiterpene emissions from Amazonian soils. The code used for including a plant functional type-dependent temperature sensitivity Nat. Commun. 9, 2226 (2018). can be found in the Supplementary Information. 28. Sarkar, C. et al. PTR-TOF-MS eddy covariance measurements of isoprene and monoterpene fluxes from an eastern Amazonian rainforest. Atmos. Chem. Phys. 20, 7179–7191 (2020). Received: 14 April 2023; Accepted: 8 December 2023; 29. Fares, S. et al. Observations of Diurnal to Weekly Variations of Monoterpene- Dominated Fluxes of Volatile Organic Compounds from Mediterranean Forests: Implications for Regional Modeling. Environ. Sci. Technol. 47, 11073–11082 (2013). 30. Saunois, M. et al. The Global Methane Budget 2000–2017. Earth Syst. Sci. Data 12, 1561–1623 (2020). References 31. Atkinson, R. Atmospheric chemistry of VOCs and NOx. Atmos. Environ. 34, 1. Loreto, F. & Schnitzler, J.-P. Abiotic stresses and induced BVOCs. Trends 2063–2101 (2000). Plant Sci. 15, 154–166 (2010). 32. Niinemets, U., Seufert, G., Steinbrecher, R. & Tenhunen, J. D. A model 2. Niinemets, Ü. Mild versus severe stress and BVOCs: thresholds, priming and coupling foliar monoterpene emissions to leaf photosynthetic characteristics consequences. Trends Plant Sci. 15, 145–153 (2010). in Mediterranean evergreen Quercus species. New Phytol. 153, 257–275 3. Loreto, F., Pinelli, P., Manes, F. & Kollist, H. Impact of ozone on monoterpene (2002). emissions and evidence for an isoprene-like antioxidant action of monoterpenes 33. Grote, R. et al. Process-based modelling of isoprenoid emissions from emitted by Quercus ilex leaves. Tree Physiol. 24, 361–367 (2004). evergreen leaves of Quercus ilex (L.). Atmos. Environ. 40, 152–165 (2006). 4. Jardine, K. J. et al. Monoterpene ‘ thermometer’ of tropical forest-atmosphere 34. Bourtsoukidis, E., Williams, J., Kesselmeier, J., Jacobi, S. & Bonn, B. From response to climate warming: Monoterpene ‘ thermometer’. Plant Cell Environ. emissions to ambient mixing ratios: online seasonal field measurements of 40, 441–452 (2017). volatile organic compounds over a Norway spruce-dominated forest in central 5. Rosenkranz, M., Chen, Y., Zhu, P. & Vlot, A. C. Volatile terpenes – mediators Germany. Atmos. Chem. Phys. 14, 6495–6510 (2014). of plant‐to‐plant communication. Plant J. 108, 617–631 (2021). 35. Tarvainen, V. et al. Temperature and light dependence of the VOC emissions 6. Peñuelas, J. & Staudt, M. BVOCs and global change. Trends Plant Sci. 15, of Scots pine. Atmos. Chem. Phys. 5, 989–998 (2005). 133–144 (2010). 36. Fares, S. et al. Biogenic emissions from Citrus species in California. Atmos. 7. Riipinen, I. et al. The contribution of organics to atmospheric nanoparticle Environ. 45, 4557–4568 (2011). growth. Nat. Geosci 5, 453–458 (2012). 37. Lim, J.-H. et al. Seasonal variations of monoterpene emissions from Pinus 8. Junninen, H. et al. Terpene emissions from boreal wetlands can initiate densiflora in East Asia. Chemosphere 73, 470–478 (2008). stronger atmospheric new particle formation than boreal forests. Commun. 38. Hakola, H., Rinne, J. & Laurila, T. The hydrocarbon emission rates of tea- Earth. Environ. 3, 93 (2022). leafed willow (Salix phylicifolia), silver birch (Betula pendula) and European 9. Kavouras, I. G., Mihalopoulos, N. & Stephanou, E. G. Formation of aspen (Populus tremula). Atmos. Environ. 32, 1825–1833 (1998). atmospheric particles from organic acids produced by forests. Nature 395, 39. Schade, G. W., Goldstein, A. H. & Lamanna, M. S. Are monoterpene 683–686 (1998). emissions influenced by humidity? Geophys. Res. Lett. 26, 2187–2190 (1999). 10. Zhang, H. et al. Monoterpenes are the largest source of summertime organic 40. Staudt, M. et al. Seasonal and diurnal patterns of monoterpene emissions from aerosol in the southeastern United States. Proc. Natl. Acad. Sci. USA 115, Pinus pinea (L.) under field conditions. Atmos. Environ. 31, 145–156 (1997). 2038–2043 (2018). 41. Laffineur, Q. et al. Impact of diffuse light on isoprene and monoterpene 11. Petäjä, T. et al. Influence of biogenic emissions from boreal forests on emissions from a mixed temperate forest. Atmos. Environ. 74, 385–392 (2013). aerosol–cloud interactions. Nat. Geosci. 15, 42–47 (2022). 42. Seco, R. et al. Ecosystem‐scale volatile organic compound fluxes during an 12. Benjamin, M. T. & Winer, A. M. Estimating the ozone-forming potential of extreme drought in a broadleaf temperate forest of the Missouri Ozarks urban trees and shrubs. Atmos. Environ. 32, 53–68 (1998). (central USA). Glob. Change Biol. 21, 3657–3674 (2015). 13. Guenther, A. B. et al. The Model of Emissions of Gases and Aerosols from 43. Filella, I., Wilkinson, M. J., Llusià, J., Hewitt, C. N. & Peñuelas, J. Volatile Nature version 2.1 (MEGAN2.1): an extended and updated framework for organic compounds emissions in Norway spruce (Picea abies) in response to modeling biogenic emissions. Geosci. Model Dev. 5, 1471–1492 (2012). temperature changes. Physiol. Plant 130, 58–66 (2007). 14. Guenther, A. et al. Estimates of global terrestrial isoprene emissions using 44. McKinney, K. A., Lee, B. H., Vasta, A., Pho, T. V. & Munger, J. W. Emissions MEGAN (Model of Emissions of Gases and Aerosols from Nature). Atmos. of isoprenoids and oxygenated biogenic volatile organic compounds from a Chem. Phys. 6, 3181–3210 (2006). New England mixed forest. Atmos. Chem. Phys. 11, 4807–4831 (2011). 15. Staudt, M. & Lhoutellier, L. Monoterpene and sesquiterpene emissions from 45. Geron, C. D., Daly, R. W., Arnts, R. R. & Guenther, A. B. & Mowry, Fred. L. Quercus coccifera exhibit interacting responses to light and temperature. Canopy level emissions of 2-methyl-3-buten-2-ol, monoterpenes, and Biogeosciences 8, 2757–2771 (2011). sesquiterpenes from an experimental Pinus taeda plantation. Sci. Total 16. Blonder, B. & Michaletz, S. T. A model for leaf temperature decoupling from Environ. 565, 730–741 (2016). air temperature. Agri. Forest Meteorol. 262, 354–360 (2018). 46. Tang, J. et al. Challenges in modelling isoprene and monoterpene emission 17. Tingey, D. T., Manning, M., Grothaus, L. C. & Burns, W. F. Influence of Light dynamics of Arctic plants: a case study from a subarctic tundra heath. and Temperature on Monoterpene Emission Rates from Slash Pine. Plant Biogeosciences 13, 6651–6667 (2016). Physiol. 65, 797–801 (1980). 47. Warneke, C. et al. Biogenic emission measurement and inventories 18. Pozzer, A. et al. Simulation of organics in the atmosphere: evaluation of determination of biogenic emissions in the eastern United States and Texas EMACv2.54 with the Mainz Organic Mechanism (MOM) coupled to the and comparison with biogenic emission inventories. J. Geophys. Res. 115, ORACLE (v1.0) submodel. Geosci. Model Dev. 15, 2673–2710 (2022). D00F18 (2010). 19. Makowski, D. et al. Systematic review of meta-analyses to assess the impacts of 48. Nagalingam, S., Seco, R., Kim, S. & Guenther, A. Heat stress strongly induces farming practices - A methodological framework. https://osf.io/byuw9. https:// monoterpene emissions in some plants with specialized terpenoid storage doi.org/10.31219/osf.io/byuw9 (2021). structures. Agri. Forest Meteorol. 333, 109400 (2023). 20. Hellén, H. et al. Sesquiterpenes dominate monoterpenes in northern wetland 49. Bäck, J. et al. Chemodiversity of a Scots pine stand and implications for emissions. Atmos. Chem. Phys. 20, 7021–7034 (2020). terpene air concentrations. Biogeosciences 9, 689–702 (2012). 21. Angot, H. et al. Biogenic volatile organic compound ambient mixing ratios 50. Stone, D., Whalley, L. K., Heard, D. E. & Tropospheric, O. H. and HO2 and emission rates in the Alaskan Arctic tundra. Biogeosciences 17, 6219–6236 radicals: field measurements and model comparisons. Chem. Soc. Rev. 41, (2020). 6348 (2012). COMMUNICATIONS EARTH & ENVIRONMENT | (2024)5:23 | https://doi.org/10.1038/s43247-023-01175-9 | www.nature.com/commsenv 9 ARTICLE COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-023-01175-9 51. Vanhatalo, A. et al. Scots Pine Stems as Dynamic Sources of Monoterpene and 79. Pringle, K. J. et al. Description and evaluation of GMXe: a new aerosol Methanol Emissions. Front. For. Glob. Change 2, 95 (2020). submodel for global simulations (v1). Geosci. Model Dev. 3, 391–412 (2010). 52. Edtbauer, A. et al. Cryptogamic organisms are a substantial source and sink 80. Fountoukis, C. & Nenes, A. ISORROPIA II: a computationally efficient for volatile organic compounds in the Amazon region. Commun. Earth thermodynamic equilibrium model for K+–Ca2+–Mg2+–NH4+–Na Environ. 2, 258 (2021). +–SO42−–NO3−–Cl−–H2O aerosols. Atmos. Chem. Phys. 7, 4639–4659 53. Mäki, M., Aalto, J., Hellén, H., Pihlatie, M. & Bäck, J. Interannual and Seasonal (2007). Dynamics of Volatile Organic Compound Fluxes From the Boreal Forest 81. Karydis, V. A., Tsimpidi, A. P., Pozzer, A., Astitha, M. & Lelieveld, J. Effects of Floor. Front. Plant Sci. 10, 191 (2019). mineral dust on global atmospheric nitrate concentrations. Atmos. Chem. 54. Byron, J. et al. Chiral monoterpenes reveal forest emission mechanisms and Phys. 16, 1491–1509 (2016). drought responses. Nature 609, 307–312 (2022). 82. Tsimpidi, A. P., Karydis, V. A., Pandis, S. N. & Lelieveld, J. Global combustion 55. Langford, B. et al. Fluxes and concentrations of volatile organic compounds sources of organic aerosols: model comparison with 84 AMS factor-analysis from a South-East Asian tropical rainforest. Atmos. Chem. Phys. 10, data sets. Atmos. Chem. Phys. 16, 8939–8962 (2016). 8391–8412 (2010). 83. Pozzer, A. et al. Distributions and regional budgets of aerosols and their 56. Yáñez-Serrano, A. M. et al. Heat stress increases the use of cytosolic pyruvate precursors simulated with the EMAC chemistry-climate model. Atmos. Chem. for isoprene biosynthesis. J. Exp. Botany 70, 5827–5838 (2019). Phys. 12, 961–987 (2012). 57. Borenstein, M., Hedges, L. V., Higgins, J. P. T. & Rothstein, H. R. Introduction 84. Pozzer, A. et al. Effects of business-as-usual anthropogenic emissions on air to Meta-Analysis. (John Wiley & Sons, Ltd, 2009). https://doi.org/10.1002/ quality. Atmos. Chem. Phys. 12, 6915–6937 (2012). 9780470743386. 85. Lelieveld, J., Evans, J. S., Fnais, M., Giannadaki, D. & Pozzer, A. The 58. Hedges, L. V., Gurevitch, J. & Curtis, P. S. The Meta-Analysis Of Response contribution of outdoor air pollution sources to premature mortality on a Ratios In Experimental Ecology. Ecology 80, 1150–1156 (1999). global scale. Nature 525, 367–371 (2015). 59. Makowski, D., Piraux, F. & Brun, F. From Experimental Network to Meta- 86. Karydis, V. A. et al. Global impact of mineral dust on cloud droplet number analysis: Methods and Applications with R for Agronomic and Environmental concentration. Atmos. Chem. Phys. 17, 5601–5621 (2017). Sciences. (Springer Netherlands, 2019). https://doi.org/10.1007/978-94-024- 1696−1. 60. Bates, D., Mächler, M., Bolker, B. & Walker, S. Fitting Linear Mixed-Effects Acknowledgements Models Using lme4. J. Stat. Soft. 67, (2015). We acknowledge the EMME-CARE project from the European Union’s Horizon 2020 61. Moreno-Martínez, Á. et al. A methodology to derive global maps of leaf traits Research and Innovation Programme (grant agreement No. 856612), as well as matching using remote sensing and climate data. Remote Sens. Environ.t 218, 69–88 co-funding by the Government of the Republic of Cyprus. We acknowledge the use of the (2018). HPC cluster Aether at the University of Bremen, financed by DFG within the scope of the 62. Camps-Valls, G. et al. Global trait maps at 1km resolution. https://doi.org/10. Excellence Initiative. We acknowledge the support of BmBf project ATTO (01LK1602B). 17871/TRY.60 (2020) 63. MAO, J. & YAN, B. Global Monthly Mean Leaf Area Index Climatology, Author contributions 1981-2015. 6.581285 MB https://doi.org/10.3334/ORNLDAAC/1653 (2019) E.B. designed the research and drafted the paper. E.B., N.D., J.P., M.V., and A.M.Y-S 64. Samek, W., Wiegand, T. & Müller, K.-R. Explainable Artificial Intelligence: collected the data. E.B., G.L., D.M., V.N.M., and A.P. analyzed the data. E.B., A.P., P.C., Understanding, Visualizing and Interpreting Deep Learning Models. J.L., J.P., J.S and J.W. performed research. All authors contributed to editing the article. arXiv:1708.08296 [cs, stat] (2017). 65. Matthaios, V. N. et al. Quantifying factors affecting contributions of roadway exhaust and non-exhaust emissions to ambient PM10–2.5 and PM2.5–0.2 Competing interests particles. Sci. Total Environ. 835, 155368 (2022). The authors declare no competing interests. 66. Matthaios, V. N. et al. Factors Influencing Classroom Exposures to Fine Particles, Black Carbon, and Nitrogen Dioxide in Inner-City Schools and Their Implications for Indoor Air Quality. Environ. Health Perspect 130, Additional information Supplementary information The online version contains supplementary material 047005 (2022). available at https://doi.org/10.1038/s43247-023-01175-9. 67. Zhou, J. et al. Random Forests and Cubist Algorithms for Predicting Shear Strengths of Rockfill Materials. Appl. Sci. 9, 1621 (2019). Correspondence and requests for materials should be addressed to Efstratios 68. Wang, Y. & Witten, I. Induction of model trees for predicting continuous Bourtsoukidis. classes. (University of Economics, Faculty of Informatics and Statistics, 1997). 69. Kuhn, M. & Johnson, K. Applied Predictive Modeling. (Springer, 2013). Peer review information Communications Earth & Environment thanks Andrea 70. Noi, P., Degener, J. & Kappas, M. Comparison of Multiple Linear Regression, Ghirardo, Thomas Holst, Tao Li, and the other, anonymous, reviewer(s) for their Cubist Regression, and Random Forest Algorithms to Estimate Daily Air contribution to the peer review of this work. Primary Handling Editors: Clare Davis. A Surface Temperature from Dynamic Combinations of MODIS LST Data. peer review file is available. Remote Sens. 9, 398 (2017). 71. Jöckel, P. et al. The atmospheric chemistry general circulation model Reprints and permission information is available at http://www.nature.com/reprints ECHAM5/MESSy1: consistent simulation of ozone from the surface to the mesosphere. Atmos. Chem. Phys. 6, 5067–5104 (2006). Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in 72. Jöckel, P. et al. Development cycle 2 of the Modular Earth Submodel System published maps and institutional affiliations. (MESSy2). Geosci. Model Dev. 3, 717–752 (2010). 73. Tost, H. et al. Global cloud and precipitation chemistry and wet deposition: tropospheric model simulations with ECHAM5/MESSy1. Atmos. Chem. Phys. 7, 2733–2757 (2007). Open Access This article is licensed under a Creative Commons 74. Ouwersloot, H. G., Pozzer, A., Steil, B., Tost, H. & Lelieveld, J. Revision of the Attribution 4.0 International License, which permits use, sharing, convective transport module CVTRANS 2.4 in the EMAC atmospheric adaptation, distribution and reproduction in any medium or format, as long as you give chemistry–climate model. Geosci. Model Dev. 8, 2435–2445 (2015). appropriate credit to the original author(s) and the source, provide a link to the Creative 75. Pozzer, A. et al. Observed and simulated global distribution and budget of Commons licence, and indicate if changes were made. The images or other third party atmospheric C2-C5 alkanes. Atmos. Chem. Phys. 10, 4403–4422 (2010). material in this article are included in the article’s Creative Commons licence, unless 76. Tsimpidi, A. P., Karydis, V. A., Pozzer, A., Pandis, S. N. & Lelieveld, J. indicated otherwise in a credit line to the material. If material is not included in the ORACLE (v1.0): module to simulate the organic aerosol composition and article’s Creative Commons licence and your intended use is not permitted by statutory evolution in the atmosphere. Geosci. Model Dev. 7, 3153–3172 (2014). regulation or exceeds the permitted use, you will need to obtain permission directly from 77. Sander, R. et al. The atmospheric chemistry box model CAABA/MECCA-3.0. the copyright holder. To view a copy of this licence, visit http://creativecommons.org/ Geosci. Model Dev. 4, 373–380 (2011). licenses/by/4.0/. 78. Lelieveld, J., Gromov, S., Pozzer, A. & Taraborrelli, D. Global tropospheric hydroxyl distribution, budget and reactivity. Atmos. Chem. Phys. 16, 12477–12493 (2016). © The Author(s) 2024 10 COMMUNICATIONS EARTH & ENVIRONMENT | (2024)5:23 | https://doi.org/10.1038/s43247-023-01175-9 | www.nature.com/commsenv