High-latitude vegetation changes will determine future plant volatile impacts on atmospheric organic aerosols

High-latitude vegetation changes will determine future


INTRODUCTION
The northern high latitudes are experiencing stronger warming than the global average and this warming is reflected in observed changes to vegetation composition, plant traits, and plant productivity 1,2 , which could profoundly alter the magnitude and composition of plant-emitted biogenic volatile organic compounds (BVOCs) 3 .Warming-induced permafrost thaw could release previously locked nutrients, abating nutrient limitations, and thereby supporting enhanced plant productivity and growth in this region 4,5 .The predicted increase in atmospheric CO 2 concentrations might, however, inhibit BVOC production, so it remains unclear how BVOC emissions might respond to the fast and combined environmental changes in the high latitudes.
Constitutive BVOC emissions from plants are dominated by terpenoids 6 , with isoprene and monoterpene contributing to ~70% and ~11% of global total BVOC emissions, respectively 7 .The production of both isoprene and terpenes is expected to increase with temperature rise 3 .As terpenoid production is coupled with plant photosynthesis, factors such as the availability of water, nutrients, light, and atmospheric CO 2 concentrations that influence plant photosynthetic rates, can indirectly impact BVOC production.Plants growing under elevated CO 2 concentrations often show lower BVOC production (clear evidence for isoprene 8,9 , but contrasting evidence for monoterpenes [10][11][12] ).However, shortterm laboratory studies have shown that the CO 2 inhibition of isoprene production decreases at higher leaf temperatures 13,14 and it is not well understood whether this pattern is sustainable in the long-term.At landscape scales, BVOC emissions vary strongly in space, depending on the plant species composition as well as their emission profiles 15 .The projected increase in atmospheric CO 2 concentrations 16 , together with rising temperatures, will likely stimulate plant productivity and thus, increase BVOC-emitting leaf biomass, while simultaneously imposing varying degrees of inhibition on BVOC production 17 and altering spatial distributions of high-latitude vegetation.
Plant-emitted BVOCs participate in a series of chemical reactions in the atmosphere.The atmospheric lifetimes of the dominant BVOCs range from minutes to hours 18 and the reactivity varies even within the same chemical group, depending on the structure of the molecule and the oxidant it reacts with.Once emitted into the air, these BVOCs, typically with the unsaturated functional groups, are susceptible to oxidation by hydroxyl radical (OH), ozone (O 3 ) and nitrate radical (NO 3 ) 19 .The major oxidation pathway of isoprene is via OH since they have similar diurnal variations.This pathway produces isoprene hydroxy peroxy radical (ISOPOO) which can either be further oxidised to form semi and extremely low-volatility organic compounds (SVOCs and ELVOCs) to contribute to secondary organic aerosol (SOA) formation 20 , or undergo H-shift isomerization to recycle OH.The oxidation of isoprene by O 3 also produces SVOCs and ELVOCs, as well as other organic peroxy radicals and small molecules 19,20 .Similarly, the SOA precursors are also formed from the oxidation of monoterpenes by OH and O 3 but with about one order of magnitude larger yields than that of isoprene 20 .As NO 3 is photolyzed rapidly, the reactions of BVOCs with NO 3 mainly occur during nighttime, which produces organic nitrates, peroxy radicals, and subsequently oxidised compounds 21 .When isoprene and monoterpenes coexist, isoprene can compete with monoterpenes to consume OH and produce relatively more volatile products than oxidation from monoterpenes.leading to suppressing SOA formation 22,23 .This suppression effect could be important to understanding synergistic feedback of changing emissions of isoprene and monoterpenes in future high-latitude climates.
In the high-latitude atmosphere, anthropogenic sources of aerosols are generally lower than in more densely inhabited regions and the Arctic is known to have low ambient concentrations of cloud condensation nuclei, CCN 24,25 .Thus, the warminginduced changes in BVOC emissions may provide stronger feedback to the high-latitude climate system, through modulation of atmospheric SOA and CCN concentrations, than elsewhere on the globe.Regional and site-specific studies have revealed that concentrations of biogenic aerosols show a clear increase with temperature [26][27][28] and have further pointed to the significance of these biogenic aerosols for climate feedback.Nevertheless, the biochemical and biophysical BVOC-mediated feedbacks on climate have been under-assessed in the high latitudes 25 .Furthermore, previous estimates of high-latitude BVOC emissions are highly uncertain 6,7 , due to the scarcity of observation-based emissions data and/or underrepresented plant variations in largescale modelling, particularly in the tundra biome 3 .Here, through integrating available BVOC measurement data and enhancing details in representing high-latitude vegetation composition in a dynamic vegetation model LPJ-GUESS, we quantitatively assess future BVOC dynamics in the Arctic and boreal regions, elucidating key processes driving the trends in BVOC emissions.Through coupling LPJ-GUESS with the atmospheric chemistry transport model, we further illustrate the contribution of plant-based BVOC emissions to our climate system through SOA-CCN-climate feedback.

Future changes in isoprene and monoterpene emissions
We explore both historical and future emissions of the dominant BVOCs, isoprene and monoterpenes, using a dynamic vegetation model, LPJ-GUESS 29 .The model, driven by climate data, simulates plant competition and vegetation composition change, as well as plant and soil biochemical processes in response to changing environmental conditions (see Methods).Observation-based BVOC emission rates and temperature response curves, together with a detailed representation of tundra plant functional types (PFTs) in LPJ-GUESS have allowed us to better represent shrub and arctic BVOC emissions 3,30 .Evaluation details of modelled plant productivity, leaf area index (LAI) and aboveground carbon over different historical periods can be found in Supplementary Figs.1-3 and the modelled ecosystem-level BVOC emissions are compared with available observational datasets across the Arctic and boreal regions Supplementary Table 1.Overall, the model can capture the spatial and temporal changes of these variables at different historical periods.
For the future period (2015-2100), we select climate scenarios from three General Circulation Models (GCMs), following five different Shared Socioeconomic Pathways (SSPs) under the CMIP6 framework 16 , and implement bias-correction of these climate predictions for temperature, precipitation, and radiation before using them to drive LPJ-GUESS (see Methods).Compared with the emission scenarios of Representative Concentration Pathways (RCPs) used in CMIP5, the SSPs in CMIP6 comprise five narratives to describe socioeconomic trends, and both how and if climate change mitigation efforts can reach the radiative forcings in the RCPs.The scenario names are a combination of the SSP narratives with the targeted RCP's radiative forcing.The resulting 15 scenarios (see the list of these scenarios in the Supplementary Table 2 and future anomalies of temperature, precipitation, and radiation in the Supplementary Fig. 4) represent five levels of (from SSP1-SSP5) varying greenhouse gas projections and three GCMs with different climate sensitivities in the study region (tundra and boreal biomes based on RESOLVE Ecoregions 2017 31 ).
We select these 3 GCMs from the full set of CMIP6 GCMs to briefly represent high, median, and low projected temperature increases over the study region.The average temperature increases by the year 2100 range from 1 to 12 °C across these 15 scenarios by the year 2100 and the warmer it gets, the stronger the increase in annual precipitation, with a maximum projected increase of 50% by the end of this century (Supplementary Fig. 4).These 15 scenarios also show both increases and decreases of incoming shortwave radiation.We use the LPJ-GUESS outputs driven by these 15 climate scenarios (hereafter "standard runs") to explore future BVOC emissions.
We simulate a clear increase in the annual, areal total emissions of isoprene (74-120% by the year 2100) and monoterpenes (11-36% by the year 2100) for the runs driven by the CanESM5 output (with the highest projected temperature increases) under different SSPs (Fig. 1d, h).For the simulations driven by scenarios with a more moderate temperature increase (such as MRI-ESM2_0 with an average temperature increase of 2.8 °C across 5 SSPs and GFDL-ESM4 with an average temperature increase of 2.4 ο C across 5 SSPs, see temperature anomalies in Supplementary Fig. 4a), total isoprene emissions show moderate increases (Fig. 1d), and total monoterpene emissions show clear decreases (see solid lines in Fig. 1h).We found that the spread in modelled total isoprene emissions across the three chosen GCMs is greater than the spread seen across the 5 SSPs from one GCM.For the monoterpenes, the spread across SSPs is more comparable to that across GCMs.
Within each grid cell, LPJ-GUESS simulates a mix of different PFTs within the canopy, and Fig. 2 shows the PFTs with the largest LAI within each grid cell, i.e., the dominant PFT.Spatially, isoprene emissions significantly increase in many regions, with the largest increasing trends simulated in regions where the dominant PFTs shift strongly (Fig. 2).The projected shifts include the replacement of boreal needle-leaved evergreen trees with broad-leaved deciduous trees (shift from PFT BNE to PFT IBS in Fig. 2) in northern Canada and western Russia, and a northward movement of boreal needle-leaved evergreen trees replacing herbaceous vegetation and shrubs in eastern Russia, Alaska, and north-eastern Canada (see PFT BNE in these regions in Fig. 2b-e).These northward migrated needle-leaved evergreen trees mainly emit monoterpenes, but can also emit a small magnitude of isoprene, and these trees are more productive than previously dominated herbaceous vegetation and shrubs, likely leading to the increase of isoprene emissions (see Supplementary Fig. 5).In the High Arctic, shrub abundance increases strongly, especially under CanESM5 SSP585 (see PFT HSS in Fig. 2 and latitudinal fractions of each PFT in Supplementary Fig. 8b).These modelled PFT shifts are in agreement with predictions based on different approaches (such as the machinelearning based ecological niche model 28 , and analysis of long-term satellite and air photos 32 ) and consistent with paleo-records of warm periods 33 (see more detailed discussion in the Section of "Model Uncertainties" in the Supplementary).
Compared with the increasing trends seen under SSP119, the modelled isoprene emissions under SSP585 in Scandinavia show decreasing trends (Fig. 1), which might be linked to the unchanged dominant PFTs, in combination with strong CO 2 inhibition of isoprene production 34 .Under SSP585, the atmospheric CO 2 concentrations reach up to 1100 ppm by the end of the 21st century, and from the lowest to the highest CO 2 increases (i.e., from SSP1 to SSP5), we can see a gradual change of the isoprene trends from a small increase to a wide decrease trends (see Fig. 1 and Supplementary Fig. 6) in the same region.For monoterpenes, the largest increasing trends occur in northern Canada and Russia (mainly for SSP585), where boreal needleleaved evergreen (in Canada) and needle-leaved deciduous (in Russia) trees with relatively high emission capacities replace the isoprene-emitting grass PFTs (GRT and C3G) in the simulations.In the southernmost study regions, we observe a clear decrease in monoterpene emissions, especially under SSP585 (Fig. 1f), as broad-leaved, isoprene-emitting, deciduous trees replace monoterpene-emitting boreal needle-leaved trees.These unfavourable vegetation shifts for monoterpenes accompany high atmospheric CO 2 increases in this climate scenario.In general, the predicted changes in isoprene and monoterpene emissions vary among climate scenarios (Fig. 1 and Supplementary Fig. 6), and show regionally varying responses linked to shifts in the dominant vegetation.These results motivate us to further assess the impacts of these individual processes (i.e., climate, atmospheric CO 2 concentration and vegetation changes).

Key processes regulating future BVOC trends
Changing temperature and atmospheric CO 2 concentrations can exert direct impacts on BVOC synthesis 9,35,36 .In addition, atmospheric CO 2 concentrations, soil nitrogen, water availability, and climate conditions not only alter photosynthetic rates, but also vegetation dynamics (including plant growth and competition, migration, and mortality) indirectly influencing BVOC emission magnitudes and composition.We group these factors that influence high-latitude BVOC emissions into four categories, i.e., climate (including temperature, precipitation and radiation), nutrient availability (mainly nitrogen, N, limitation in the high latitudes), vegetation dynamics, and atmospheric CO 2 concentration.As some factor depends on others, such as, climate, CO 2 and nutrients can all influence vegetation dynamics, we separately checked how CO 2 , nutrient and climate influence BVOCs emissions and also together with these impacts through vegetation dynamics in regional BVOC emission changes, using a set of factorial simulations (Table 1) based on climate forcing data from GFDL-ESM4 and CanESM5's outputs with the lowest and highest emission scenarios (i.e., SSP119 and SSP585, respectively).The modelled temperature changes from these selected climate scenarios (2.6 and 6.4 °C temperature increases averaged over 2015-2100, respectively) are representative of the range of temperature changes across the 15 scenarios.We design four factorial experiments, in each of which we change one process at a time: (1) Ecosystem dynamics unaffected by CO 2 increase, represented with a constant CO 2 concentration (using the value for 2014) for the future period (hereafter, termed "noCO2"); (2) BVOC dynamics unaffected by further CO 2 inhibition, represented with constant CO 2 inhibition impacts on BVOC production (using the inhibition level for 2014) for the future period (noCO2Inhibition); (3) A removal of N limitation, achieved by adding 50 kg N ha −1 yr −1 to the annual nitrogen deposition input fields in the future period (noNlim).This N addition corresponds to what has been implemented in forest N fertilisation trials in Davies-Barnard, et al. 37 ; and (4) Ecosystem dynamics unaffected by climate change, achieved by using the monthly averages of climate drivers from the period 2005-2014 for driving ecosystem processes in the future period (noVegDym).In this simulation, the predicted temperature changes still affect BVOC production, but we remove future climate impacts on vegetation dynamics.Subsequently, we calculate the differences between the standard and factorial simulations (see Table 1) to tease apart the relative importance of CO 2 fertilisation, CO 2 inhibition on BVOC production, N limitation, and vegetation changes as determinants of spatial and temporal patterns of future BVOC emissions (Fig. 3).
Under the low CO 2 emission scenario (CanESM5 SSP119), the overall positive trend in isoprene emissions is largely driven by vegetation changes (Fig. 3a, e).The small increasing emission trends by CO 2 inhibition in both isoprene and monoterpene emissions are driven by decreasing atmospheric CO 2  6, and the probability density of significant trends for all five SSPs is shown in Supplementary Fig. 7.
concentrations towards the end of the century in this scenario (see CO 2 inhibition in Fig. 3c, m).Overall, the impacts from CO 2 fertilisation and N limitation are limited for both isoprene and monoterpenes under CanESM5 SSP119 (Fig. 3b, d, l, n).
Under the high CO 2 emission scenario (CanESM5 SSP585) with associated stronger warming and large increases in N deposition, the positive trends in isoprene emissions are associated with CO 2 fertilisation of photosynthesis and vegetation changes, although high CO 2 concentration simultaneously inhibits BVOC production (Fig. 3g, h, j).Climate warming-induced vegetation changes promote the overall positive trend in isoprene (Fig. 3j), but not in monoterpene emissions, as seen in the simulation of positive and negative impacts under CanESM5 SSP585 (Fig. 3t).The impacts from N limitation are again, very small, and likely linked to the increased N deposition during this century in CanESM5 SSP585 and increased sources of mineral N from warming soil.
We quantified the potential synergistic effects of four processes under different scenarios.Our findings suggest that the synergistic effects on BVOC emissions are typically less significant than the effects of each individual process.We also conduct the same factorial experiments under climate scenarios GFDL-ESM4 SSP119 and GFDL-ESM4 SSP585, and we see similar patterns in terms of important controlling processes (see Supplementary Fig. 9) for future isoprene and monoterpene emissions.
BVOC impacts on the regional atmosphere The modelled high-latitude BVOC emissions, together with vegetation status (including LAI, and vegetation cover fraction) of the years 2009 and 2100 from standard runs driven by CanESM5 SSP119 and CanESM5 SSP585 are fed into the global chemistry transport model, version 5 (TM5, Bergman et al. 38 ).The modelled emissions from these two selected years, i.e., 2009 and 2100, are representative for the decadal average monthly emissions (see Supplementary Fig. 10) and we, therefore, use the outputs from these two years to represent historical and future periods.Then, to provide all the input data needed at a global scale and to ensure that the changes to the atmospheric chemistry originate from high-latitude changes only, the emission data, LAI, and vegetation cover south of the study domain were all set to values from the same standard LPJ-GUESS global run from the year 2009.Furthermore, for the year 2100, we use LPJ-GUESS outputs from standard runs of CanESM5 SSP585 and CanESM5 SSP119, as well as the outputs from different factorial simulations as inputs to TM5.In the result below, we present the TM5 outputs from two of the factorial experiments: i.e., noVegDym and noCO2Inhibition, where their associated processes contribute most to the simulated BVOC trends based on the averaged values and numbers of grid cells with significant trends (see the 3 rd and 5 th columns in Fig. 3, the trends from the processes of CO 2 inhibition and vegetation changes) driven by CanESM5 SSP585 and CanESM5 SSP119 (see Methods for the detailed setup for TM5).TM5 is used to quantify the impacts of isoprene and monoterpene emissions on surface SOA concentrations (SOA surf ), SOA optical depths at 550 nm (OD550 soa ), and aerosol optical depths at 550 nm (OD550 aer ).The SOA formation in TM5 is mainly in proportion to these two dominant BVOC emissions.Then, the aerosol optical depths describe how much the light penetration through the atmosphere is prevented by aerosols.550 nm is at or near the energy peak of the solar radiation spectrum, which makes it a good representative of the incoming shortwave radiation spectrum, as applied in previous studies (e.g., refs. 39,40).OD550 soa is the optical depth considering only the SOA component, while OD550 aer is the optical depth considering all the aerosol components including SOA, sulfate, methanesulfonic acid, ammonium nitrate, black carbon, primary organic aerosols, sea salt, and mineral dust.
Our results demonstrate the important role of climate changedriven vegetation changes in regulating the spatial patterns of BVOC impacts on regional atmospheric aerosols.The increased BVOC emissions have largely contributed to the increase in SOA surf for a major part of the study region, with a smaller area of increase for CanESM5 SSP119 than for CanESM5 SSP585 (Figs. 4g and 5g).Notably, the direct impact of temperature on SOA yields is not yet accounted for in the current model version.Nevertheless, the temperature impact on the chemical reactions has been considered as an indirect impact of temperature on SOA formation.Under both scenarios, we see up to a 2.7-fold increase of SOA surf in northern Canada and Russia in the standard run.The simulation without vegetation responses to climate change (noVegDym) results in a considerably lower increase in SOA surf (Figs.4h and 5h).The standard runs with vegetation changes show stronger increases over a larger area in aerosol optical depth (OD550 soa and OD550 aer : 41% and 4.9% increase under CanESM5 SSP119, Fig. 4j, m, and 29% and 4.1% increase under CanESM5 SSP585, Fig. 5 j,  m, respectively).Without vegetation changes (noVegDym), only a limited increase of OD550 soa and OD550 aer in a small region in central and eastern Canada is simulated.As the changes in OD550 soa and OD550 aer affect the reflection and absorption of visible light in the atmosphere, an increase in these values means that a reduced amount of radiation is received at the surface.
The strong spatial linkages that we show for vegetation shifts, BVOC dynamics, and SOA/aerosol changes are not captured when using static vegetation distributions for future conditions 6 .Under the warmest scenario (i.e., CanESM5 SSP585), the widespread replacement of boreal needle-leaved trees with broad-leaved deciduous trees (Fig. 2e) in southern Finland and western Russia clearly contributes to the local reductions of SOA surf and optical depths (Fig. 5g, j, m).When future CO 2 inhibition of BVOC production is excluded under CanESM5 SSP585, TM5 estimates up to 10-fold increases in SOA surf (mainly in the high latitudes), up to a 1.2-fold increase in OD550 soa , and up to an 18% increase in OD550 aer in the year 2100, as compared to 2009.Similar spatial patterns with slightly smaller magnitudes are observed when CO 2 inhibition is excluded only for monoterpenes (see Supplementary Figs.11 and 12).Under SSP119, the effects of atmospheric CO 2 on regional SOA and aerosols via BVOCs are limited.Our results emphasise the importance of monoterpenes in influencing SOA and aerosol yields (as also pointed out in previous studies 38 ), and show that we need to understand whether CO 2 inhibition affects monoterpene production in high-latitude plants similarly to those plants that were studied in the literature 9,35,41 .
The simulated changes in BVOCs and SOA lead to increases in CCN concentrations at a supersaturation of 1.0% near the surface, mainly for the Arctic region (see Figs. 4p-r and 5p-r), which indicates the potential for enhanced formation of low-level clouds.CCN (1.0%) represents the number concentration of particles larger than 50 nm in diameter, so it is sensitive to new particle formation and growth, which are affected by the gas precursors, ELVOCs and SVOCs 20 .There is a northward shift of increased CCN (1.0%), suggesting the potential for more clouds at the high latitudes and fewer at mid-latitudes.Without CO 2 inhibition under CanESM5 SSP585, CCN over Greenland and eastern Canada increased by 16% (Fig. 5r), which highlights the strong expected impacts of CO 2 inhibition on BVOC emissions.
Based on the modelled changes in SOA and aerosol optical depth (AOD) between the years 2100 and 2009, we further Table 1.Overview of factorial runs conducted with LPJ-GUESS following CanESM5 and GFDL-ESM4 SSP119 and SSP585.2).Based on parameters from three different references 27,40,42 , we find that the radiative forcing estimates vary greatly depending on the method used.Yli-Juuti et al. 27 derived the relationships between temperature, aerosol loading, and radiative forcing feedback based on summertime (July-August) field measurements in a boreal ecosystem, and based on this study, the differences in the estimated areal averaged radiative forcing between standard and noVegDym runs are smaller than the estimations from the other two studies 40,42 with parameterisations derived from annual data at the global scale.The overall cooling feedbacks are strongest for the noCO2Inhibition run following the CanESM5 SSP585 scenario.At the grid cell level, the strongest increase and decrease in radiative forcing based on Yli-Juuti et al. 27 in the standard run are -2.09and 0.79 W m −2 , respectively, and these two values are much weaker in the noVegDym run (which are -1.52 and 0.54 W m −2 , correspondingly), highlighting the importance of warming-induced vegetation changes on regulating local climatic impacts.Under CanESM5 SSP119, the changes in areal averaged radiative forcing over July and August are -0.44,-0.34 and -0.43 W m −2 for standard, noVegDym, and noCO2Inhibition simulations, respectively.The strongest cooling feedback can be seen with warming-induced vegetation changes included in the standard run under this climate scenario.

DISCUSSION
Our results illustrate that vegetation changes in a warmer climate play a crucial role in shaping future BVOC feedback to Fig. 3 Histograms of trends in the modelled isoprene (first row for CanESM 5 SSP119 and second row for CanESM 5 SSP 5 8 5 ) and monoterpene (third row for CanESM 5 SSP119 and fourth row for CanESM5 SSP 5 8 5 ) emissions.From standard runs (a, f, k, p) and from four investigated processes (b, g, l, q CO 2 fertilization; c, h, m, r CO 2 inhibition; d, i, n, s N limitation; and e, j, o, t Vegetation changes) over the period 2015-2100.x-axis shows the number of grid cells and y-axis shows the values of estimated trends.The values used in the trend analysis of these four investigated processes are based on the differences of every two runs (see how the impacts from each process are extracted in Table 1).Trends are analysed using the Mann-Kendall test, and the trends for all grid cells are marked as grey, while both significant positive and significant negative trends (p < 0.05) are indicated with red and blue bars, respectively.Sig.Significant.
Fig. 4 The inputs to and outputs from TM 5 using CanESM 5 SSP 119 .The first and second rows show LPJ-GUESS simulated isoprene (ISO) and monoterpene (MT) emission changes between 2100 and 2009.The emissions from the year 2100 are driven by CanESM5 SSP119.The third to the sixth rows show the TM5 simulated ratio in changes to surface SOA concentration (SOA surf ); optical depth of SOA at 550 nm (OD550 soa ); optical depth of aerosol at 550 nm (OD550 aer ), and CCN concentrations at a supersaturation of 1.0% near the surface (CCN1.0).From left to right, we show the TM5 results fed with BVOC inputs from three LPJ-GUESS runs: the standard run (the 1st column, a, d changes in ISO and MT emission, g, j, m, p changes in SOA surf , OD550 soa , OD550 aer and CCN1.0), the noVegDym run (the 2nd column, b, e changes in ISO and MT emission, h, k, n, q changes in SOA surf , OD550 soa , OD550 aer and CCN1.0) and the noCO2Inhibition run (the 3rd column, c, f changes in ISO&MT emission, i, l, o, r changes in SOA surf , OD550 soa , OD550 aer and CCN1.0).The colour bars used for these three columns are kept the same for each corresponding output.
Fig. 5 The inputs to and outputs from TM 5 using CanESM5 SSP 5 8 5 .The first and second rows show LPJ-GUESS simulated isoprene (ISO) and monoterpene (MT) emission changes between 2100 and 2009.The emissions from the year 2100 are driven by CanESM5 SSP585.The third to the sixth rows show the TM5 simulated ratio in changes to surface SOA concentration (SOA surf ); optical depth of SOA at 550 nm (OD550 soa ); and optical depth of aerosol at 550 nm (OD550 aer ), and CCN concentrations at a supersaturation of 1.0% near the surface (CCN1.0).From left to right, we show the TM5 results fed with BVOC inputs from three LPJ-GUESS runs, which are the standard run (the 1st column, a, d changes in ISO and MT emission, g, j, m, p changes in SOA surf , OD550 soa , OD550 aer and CCN1.0), noVegDym run (the 2nd column, b, e changes in ISO and MT emission, h, k, n, q changes in SOA surf , OD550 soa , OD550 aer and CCN1.0) and noCO2Inhibition run (the 3rd column, c, f: changes in ISO&MT emission, i, l, o, r changes in SOA surf , OD550 soa , OD550 aer and CCN1.0).The colour bars used for standard and noVegDym runs are kept the same for each corresponding output.atmospheric chemistry and climate.Most previous assessments of feedback between the land surface and atmosphere in the high latitudes, that involve vegetation changes, have focused on changes in surface albedo 28,43 and increases in atmospheric water vapour 44 , but our study clearly demonstrates the strong and regionally diverse feedback of vegetation changes on our climate through the BVOC-SOA pathway.
The model we used in this study has been recently developed by including a range of high-latitude processes: such as a variety of shrubs, permafrost, wetland, nitrogen cycle, and by parameterising high-latitude BVOC emissions 3,30,45 .These improvements make it possible to assess the high-latitude BVOC emissions under different future scenarios.The warming-induced widespread increase of broad-leaved deciduous trees at the expense of boreal evergreen needle-leaved trees in the boreal region suppressed the emissions of monoterpenes and thereby, SOA formation, causing a regional warming feedback (up to 0.79 W m −2 for summer months in the standard run, driven by CanESM5 SSP585).In the High Arctic, the increased abundance of shrubs and the northward advance of boreal needle-leaved trees in northern Canada and Siberia contributed to an increase of surface SOA, resulting in an increase of up to 41% in SOA optical depth and likely leading to a cooling feedback (with the strongest local cooling of -2.25 and -2.09W m −2 for standard runs driven by CanESM5 SSP119 and SSP585, respectively).Currently, the Arctic features a lack of aerosol particles for cloud formation 24 and the northward shifts of vegetation bring in a new, important aerosol source: plant-emitted BVOCs.This 'new' source of aerosols might enhance cloud formation in this region, as we can see an overall increase of CCN under different runs, and can particularly increase the lowlevel cloud formation in the Arctic 25 .During the growing season, the enhanced cloud cover leads to cooling feedbacks through scattering shortwave radiation that counteracts warming feedbacks associated with the re-emission of longwave radiation received from the land surface 24,25 .Furthermore, the northward shifts of woody plants as well as changes to PFT compositions simulated by LPJ-GUESS for future periods might be overestimated as constraints such as seed dispersal are not considered in the model.In spite of this, our modelled vegetation responses are consistent with published experimental and other modelling studies (see more detailed discussion in the Supplementary Model Uncertainties section).Our estimations showed the overall cooling feedback from warming-induced changes in BVOC and SOA (negative aerosol radiative forcing, -0.44 and -0.45 W m −2 for CanESM5 SSP119 and CanESM5 SSP585, respectively, when comparing the year 2100 and the year 2009) during summer months, and these regional radiative forcing values are comparable to the effective radiative forcing (ERF) caused by aerosolradiation interactions over the industrial era globally (-0.3 W m −2 for 1750-2014) 46 .At an annual scale, the estimated aerosol regional RF (-0.07 and -0.08 W m −2 for CanESM5 SSP585 and CanESM5 SSP119, respectively, based on the averaged RF calculated from refs. 40and 42 ) are of a similar magnitude to the impacts from the stratospheric water vapour (0.05±0.05 W m −2 46 ).We note that our simplistic methods of estimating aerosol impacts on radiative forcing are strongly linked to the assumptions originating from the references used 27,40,42 .The overall net feedback from the vegetation-BVOC-SOA pathway thus needs to be comprehensively evaluated using a more comprehensive Earth System Model to fully consider the negative feedback from scattering and cloud formation from SOA, and also the positive feedbacks from increased longwave radiation emitted from lowlevel clouds 25,47 .Additionally, the spread of climate sensitivity from different GCMs needs to be reduced to better understand this pathway of feedback.
Under high CO 2 emission scenarios, the climate warminginduced increase of natural aerosols from plant BVOC emissions is largely constrained by CO 2 inhibition of BVOC (mainly monoterpene) production, which means that future anthropogenic CO 2 increase might provide an indirect positive feedback to the climate through this inhibition.Meanwhile, constrained BVOC emissions under future high atmospheric CO 2 concentration might increase hydroxyl radical (OH) concentrations, and then reduce the lifetime of methane, providing an indirect negative feedback to the climate 48 .The potential for strong atmospheric feedback associated with plant BVOC responses to the increasing atmospheric CO 2 concentrations makes the need for more leafand ecosystem-level observations to unveil the mechanisms for the decoupling between photosynthesis and BVOC production under elevated CO 2 .This is especially relevant in the context of climate change as the CO 2 inhibition of isoprene emission may be suppressed at higher temperature 13,14 , adding uncertainty to the long-term impacts of rising CO 2 concentrations on isoprene emissions.Furthermore, it is still unclear whether elevated CO 2 inhibits monoterpene production to the same degree as isoprene production.The empirical function used in LPJ-GUESS to assess the CO 2 impact on terpenoid production 34 was derived from a limited number of observational studies on trees.Whether the BVOC emissions of low-statured Arctic plants respond similarly to CO 2 is unknown.We currently have no published data that enable a quantification of the low-statured plants' BVOC responses to the surrounding CO 2 .
As the current TM5 simulation settings neither consider the impacts from future changes in meteorology, nor other surface emissions except isoprene, monoterpenes emissions and vegetation cover, the current TM5 setup allows us to single out the 'isolated' impacts from plant-emitted isoprene and monoterpenes alone, regardless of interactions with future changes in other factors (see "Model uncertainties" in the Supplementary).Laboratory experiments have shown that higher temperatures can The values in the brackets are minimum and maximum radiative forcing in 2100 compared to the year 2009 over the study region.Notably, the parameters based on the ref. 27is mainly based on the study of summertime field measurement in a boreal ecosystem, and the parameters are based on refs. 40,42are derived from annual data at the global scale.
decrease SOA yield 49,50 , as under warmer conditions, BVOC reaction products remain volatile for a longer time.Furthermore, isoprene in the gas mixture, through OH scavenging can suppress new particle formation 51 and SOA formation 22 .These processes are not yet included in the large-scale model TM5, so we can only speculate that the scenarios with the largest predicted temperature increases may overestimate the SOA yield.Globally, we expect a reduction in anthropogenic emissions, such as such as SO 2 and NOx as presented in the (See IPCC Technical Summary report, Fig. TS.4 in ref. 52 ) under different SSPs.These changes mainly occur south of the study domain and can influence aerosol seeds and chemistry in the Arctic.It might be the case that the decrease of SO 2 dampens aerosol formation as well as potential cooling impacts in the high-latitude regions, but we are not sure about the magnitude yet, as this is influenced by the background concentration (clean environment) in the high-latitude regions as well as depending on the reduced amount of the transported SO 2 .Climate warming will also likely change the baseline emissions from open oceans, and potentially increase marine biogenic trace gas, such as dimethylsulfide (DMS), emissions, which can contribute to new particle formation 25 .With the predicted further reduction of the sea-ice extent, there will likely be more shipping and gas extraction in the Arctic, which could bring in new emissions of SO 2 , NO X and black carbon etc, which might counteract the emission reductions of these compounds further south 53 .Future studies should focus on quantification of the synergistic effects of future plant emissions from high latitudes with other anthropogenic and primary aerosol sources under dynamic changing climate conditions in a coupled Earth System Model (such as ref. 54).In summary, our results show a potentially significant feedback mechanism linking climate change-induced vegetation changes, BVOC dynamics and aerosols in the high latitudes.The negative feedback mechanism between the biosphere, aerosols, and climate concluded from observational data 55 cannot be extrapolated into the future without considering climate changeinduced vegetation changes and their impacts on emitted compounds.Our study confirms the importance of BVOCs for future atmospheric SOA concentration, optical depths and associated radiative forcing in the high latitudes.Furthermore, the study also reveals that the overall impacts largely depend on how atmospheric CO 2 concentration influences monoterpene production.The net radiative feedbacks from BVOCs need a comprehensive evaluation in order to assess (1) vegetation change-induced shifts in emission profiles and (2) the balance between aerosol shortwave cooling and aerosol longwave warming feedbacks in the high-latitude environment.

Dynamic vegetation model, LPJ-GUESS
We used the latest version of LPJ-GUESS v4.1 with the relevant developments of wetland biogeochemistry and soil physics following 56,57 .LPJ-GUESS is a dynamic ecosystem model, which simulates vegetation growth, mortality, and competition, as well as soil biogeochemistry 29 .The model has been widely used to assess water, nitrogen and carbon fluxes, as well as vegetation dynamics at regional and global scales.Plants are represented as plant functional types (PFTs) with a set of predefined bioclimatic, physiological, life history and phenological parameters that characterise specific plant growing requirements, structure and spatial distribution.For simulations in high latitudes, different levels of shrubs (high, low and prostrate), lichen, moss and wetland PFTs are specified 30 .Compared with the version of LPJ-GUESS described in Smith et al. 29 , the model version used here includes improved soil temperature calculations, allowing a better representation of soil thawing and freezing processes and influencing water availability to plants.The model also includes wetland biogeochemistry and wetland hydrology 58 .
Daily leaf-level photosynthesis is based on the simplified version of the Farquhar biochemical model 59,60 , which simulates a gradual transition between an electron-transport-limited and a rubisco-limited carbon assimilation.In the model, a fraction of the photosynthetic electron-flux is used to produce isoprene 34 and monoterpenes 61 , and the production rates of isoprene and monoterpenes are also influenced by PFT-specific emission factors, leaf temperature, seasonality and atmospheric CO 2 concentration.The PFT-specific emission factors (standardised emission rates at photosynthetically active radiation levels of 1000 µmol m −2 s −1 and at a reference temperature of 20 °C for Arctic PFTs) in the tundra region are based on the branch-or leaflevel measurement data (see details in refs. 3,30), and a stronger temperature sensitivity derived in Tang et al. 30 has been applied for Arctic PFTs.For boreal and temperate PFTs, the global temperature response curve has been used and the emission factors are based on the reference temperature of 30 °C.The model dynamically compares plant water demand with liquid water supply from the model's 15, 10-cm thick soil layers to determine transpiration rate.If the water supply cannot meet the demand, plant photosynthetic rates are downregulated 29 , which impacts the electron-flux used for isoprene and monoterpene production.In this way, the model considers substrate availability for BVOC production under drought and in the presence of frozen soil layers, but the model does not yet consider any potential effects of drought stress otherwise.

Standard runs
For simulation of the historical period in this study, we used temperature, precipitation and shortwave radiation fields from the monthly CRU-NCEP climate dataset (https://rda.ucar.edu/datasets/ds314.3/)at 0.5 ο by 0.5 ο resolution as input to the model for the period of 1901-2014.Detrended CRU-NCEP fields from the 1901-1930 period were repeated and looped through as input to the model for the period 1850-1900.For the future period (i.e., 2015-2100), we selected the climate outputs from three general circulation models (GCMs), each following five different Shared Socioeconomic Pathways (SSPs) from the latest CMIP6 project 16 to represent a range of predicted future climates (see LPJ-GUESS Standard runs in Fig. 6).The three GCMs were selected based on our analysis of temperature changes in the high latitudes, and also on the completeness of the published outputs from these five SSPs (as of September 2020) to represent a wide range of future climate changes in the study region.The monthly climate data from each GCM and SSP for the future period were bias-corrected.The biases were calculated as the difference between the monthly climate data (using delta change methods for temperature, and relative anomalies for precipitation and radiation) of the period 1985-2014 from CRU-NCEP and climate outputs from each GCM in the same period, and these biases were then only implemented to the future climate simulated by the GCM.A detailed description of the bias-correction approach can be found in Ahlström et al. 62 .The predicted anomalies in future temperature, precipitation, and radiation for the 15 scenarios (3 GCMs x 5 SSPs) are presented in Supplementary Fig. 4. All 15 standard runs share a common climate development during the historical period and start to diverge from 2015 onwards following the different future scenarios.

Factorial simulations in LPJ-GUESS
Four factorial experiments were designed to separate the impacts on BVOC dynamics from all relevant processes in the model (see Fig. 6).We have described how each factorial experiment was implemented in the main text (section: Key processes regulating future BVOC trends).Briefly, we changed one process at a time and compared the differences between the standard and factorial simulations over the period 2015-2100.The factorial simulations were conducted under the climate scenarios of CanESM5 SSP119 & SSP585 (see the results in Fig. 3), as well as GFDL-ESM4 SSP119&SSP585 (see the results in Supplementary Fig. 9).We also conducted one additional factorial experiment where only monoterpene emissions are unaffected by CO 2 inhibition (namely noMT_CO2Inhibition) under CanESM5 SSP119 and SSP585 (see Supplementary Figs.11 and 12).

Global chemistry transport model, TM5
To further assess the impacts of plant-emitted BVOCs on atmospheric aerosols and cloud condensation nuclei (CCN), the modelled leaf area index and vegetation coverage, together with isoprene and monoterpene emissions were fed into a global chemistry transport model, TM5-MP 63 .TM5-MP is a branch of TM5 with a massively parallel functionality and is now maintained by KNMI (Royal Netherlands Meteorological Institute).For simplicity, we refer to the model as TM5 throughout the text.The meteorological and surface fields driving the model were derived from ERA-Interim reanalysis datasets provided by ECMWF (European Centre for Medium-range Weather Forecasts) 64 , which are the default forcing datasets for TM5.The chemistry scheme used in this study is a modified version of CB05 (carbon bond mechanism; Yarwood et al. 65 ) with more details described in Williams et al. 63 .Aerosol processes are calculated with the modal two-moment model M7 66 .It includes seven log-normally distributed modes comprising four water-soluble modes (nucleation, Aitken, accumulation and coarse) and three insoluble modes (Aitken, accumulation and coarse).The dry diameter range of each mode is <10 nm for nucleation mode, 10 nm to 100 nm for Aitken mode, 100 nm to 1000 nm for accumulation mode, and >1000 nm for coarse mode.The CCN concentration at a given supersaturation is calculated with volume weighted κ (hygroscopicity parameter) values which depend on the aerosol composition, following the activation scheme in ref. 67 .
Originally in TM5, inputs of monthly mean natural emissions of isoprene and monoterpenes are derived from MEGANv2.1 6,7 .Then a diurnal cycle is applied to the monthly mean values.However, in this study, we substituted these monthly mean emission data with the monthly emission outputs from individual LPJ-GUESS simulation runs.The emissions of isoprene and monoterpenes from biomass burning were applied from the default inventory provided by van Marle et al. 68 without diurnal variations.Furthermore, the oceanic dimethylsulfide (DMS) emissions, the mineral dust emissions and the sea salt emissions are calculated within TM5 69 .The other natural emissions such as CO, nonmethane VOCs, NOx (NO + NO 2 ), NH 3 and SO 2 were prescribed as in van Noije et al. 70 .The anthropogenic and biomass burning emissions of gases and particles were derived for present-day conditions from the CMIP6 input4MIPs inventory 68,71 .
Once emitted, isoprene and monoterpenes can react with hydroxyl radical (OH) and ozone (O 3 ) to produce ELVOCs (extremely low volatile organic compounds) and SVOCs (semivolatile organic compounds), which can condense on particles to increase SOA mass.In addition, ELVOCs can participate in new particle formation together with sulfuric acid.These processes were recently implemented in TM5, which led to improvements in comparisons with observed aerosol concentration and satellite estimates of AOD (see Bergman et al. 38 for more details).
In this study, a horizontal resolution of 3 degrees in longitude and 2 degrees in latitude was applied in TM5.In the vertical direction, 34 hybrid-sigma levels were used.The time step was one hour.All the simulations were run for the year 2009 with a spin-up period of one year.The meteorological and surface fields except for the vegetation cover in 2009 were applied for all the simulation cases, which omitted the meteorological impacts in future scenarios when compared to the present case (see TM5 in Fig. 6).In the current version of TM5-MP, the change in vegetation cover does not directly affect the aerosol dry deposition velocity.However, it will affect the dry deposition of gas-phase species, subsequently affecting the aerosol concentration.For example, the concentration changes of O3 and SO2 can affect the SOA formation via chemical reactions related to isoprene and monoterpenes.Furthermore, the non-BVOC emission datasets applied in all simulations were from the year 2009 except those derived from LPJ-GUESS output as mentioned above.TM5 was installed and configured in CSC (Finnish IT Center for Science) Puhti.

Estimation of radiative forcings caused by SOA changes
The radiative forcings due to SOA changes in six simulation cases (including standard, noVegDym, and noCO2Inhibition under two climate scenarios CanESM5 SSP585 and SSP119) were calculated according to the methods and results described in Zhu et al. 42 , Yli-Juuti et al. 27 and Bellouin et al. 40 , respectively (Table 2).The model results in Zhu, et al. 42 showed that SOA burden can increase by 6.8% in 2100 compared to 2000 under the RCP8.5 scenario considering the combined effects of the future changes of climate, anthropogenic emissions and land use.From their results, we can derive that the changing rate of radiative forcing due to SOA burden change is -0.441W m −2 mg −1 m 2 .The radiative forcing estimated from our simulation cases was then calculated by multiplying this value by the burden change over the study domain (including the Arctic and boreal ecosystems).Yli-Juuti et al. 27 analysed multiple datasets at a boreal measurement station Hyytiälä in Finland.Using a multivariate mixed-effect model, they found that in July and August, the organic aerosol (OA) mass loading changed with temperature at a rate of 0.24 (ranging from 0.22 to 0.25) μg m −3 K −1 while the all-sky total radiative forcing feedback based on sun photometer observations is -0.47 W m −2 K −1 .We derived the radiative forcing due to SOA surface concentration change by dividing the radiative forcing feedback by OA mass loading change, which is equal to -1.96 W m −2 μg −1 m 3 .Based on this, we estimated the radiative forcing feedback for our simulation cases using the SOA surface concentration in July and August over the study domain.Bellouin et al. 40 calculated the present-day effective radiative forcing compared to pre-industrial with a theoretical formula which is a function of aerosol optical depth change.This formula with the suggested parameters was applied to calculate the radiative forcing of our simulation cases both at the global scale and for the study region.All the parameters in Bellouin et al. 40 were estimated with a confidence interval, and we only used their median values to show an estimation of radiative forcing due to SOA change.

Fig. 1
Fig. 1 LPJ-GUESS modelled trends of isoprene and monoterpene emissions for the period of 2001-2100.a, b Modelled isoprene (ISO) trends for SSP119 and SSP585, respectively.The trends are analysed based on the averaged emissions over 3 General Circulation Models (GCMs) and only significant trends (Mann-Kendall trend test, p < 0.05) are shown; c Probability density of significant trends (shown in panels a and b) in isoprene emissions with kernel smoothing distribution fit; d Time series of areal total (study region: tundra and boreal biomes) isoprene emissions for all standard runs driven by 3 GCMs following 5 SSPs; e, f Modelled monoterpene (MT) trends for SSP119 and SSP585, respectively.The trends are analysed based on the averaged emissions over 3 GCMs and only significant trends (Mann-Kendall trend test, p < 0.05) are shown; g Probability density of significant trends (shown in panels e and f) in monoterpene emissions with kernel smoothing distribution fit; h Time series of areal total monoterpene emissions for all standard runs driven by 3 different GCMs following 5 SSPs.SSP119: Shared Socioeconomic Pathway 1 reaching radiative forcing of 1.9 W/m 2 in 2100.SSP585: Shared Socioeconomic Pathway 5 reaching radiative forcing of 8.5 W/m 2 in 2100; d, h For each GCM, the solid line represents the average across 5 SSPs and the shaded area represents the predicted ranges from minimum to maximum.LPJ-GUESS runs at 0.5 ο by 0.5 ο resolution.The same trend maps made for SSP126, SSP245 and SSP370 are shown in Supplementary Fig. 6, and the probability density of significant trends for all five SSPs is shown in Supplementary Fig.7.

Fig. 2
Fig. 2 LPJ-GUESS simulated historical and future distribution of the dominant plant function types (PFTs).Distribution of the dominant plant functional types (PFTs) over the period 1971-2000 (a) and 2071-2100 (b-e) based on the modelled leaf area index.The outputs from the scenario SSP119 are shown in b, c, and the outputs from SSP585 are shown in d, e.The outputs from two GCMs are plotted separately: GFDL-ESM4 (b, d) and CanESM5 (c, e).
effects on aerosol radiative forcing among the standard, noVegDym, and noCO2Inhibition runs (Table

Fig. 6
Fig.6Model setup.The modelled monthly isoprene (ISO) and monoterpene (MT) emissions, as well as vegetation status (Veg) from LPJ-GUESS standard and factorial runs, were used as inputs for TM5.The years 2009 and 2100 were selected to represent historical and future periods in the global model TM5.ECMWF European Centre for Medium-Range Weather Forecasts.SSPs Shared Socioeconomic Pathways (SSPs).

Table 2 .
Estimated radiative forcing (W m −2 ) in 2100 in comparison to the year 2009, averaged over the study domain (i.e., the Arctic and boreal ecosystems).