Toluene early morning peaks as potential biogenic stress signals of Mediterranean forests

Abstract

Mediterranean forests emit biogenic volatile organic compounds (BVOCs) that significantly influence atmospheric chemistry and ecosystem functioning. Due to high solar radiation and the important anthropogenic emission of atmospheric pollutants in the Western Mediterranean Basin, photochemical activity is enhanced, favouring the formation of ozone and secondary organic aerosols. We measured ambient mixing ratios of volatile organic compounds (VOCs) using Proton Transfer Reaction Mass Spectrometry (PTR-MS) at Montseny Natural Park, a Mediterranean forest 60 km from Barcelona, during summers 2021–2023. Surprisingly, during heat wave conditions, a toluene early morning peak of 0.23 (interquartile range 0.18) ppbv was observed. These early morning peaks were associated with prolonged high temperatures (approximately 5 °C higher as compared to non-peak days) and persistently higher vapour pressure deficit (approximately 877 Pa higher as compared to non-peak days), suggesting a potential link to vegetation drought stress. Cross-validation of measurements with Gas Chromatography-Mass Spectrometry (GC-MS) enhanced data reliability. Analysis of the meteorology and the variability of other atmospheric pollutants at the site, allowed to attribute these peaks to a local biogenic origin. These findings suggest that early morning toluene peaks could serve as an indicator of Mediterranean forest vulnerability to climate-induced stressors, with potentially more frequent and intense peaks occurring in the future.

---

Environmental significance

In this work, we present high-resolution measurements of volatile organic compounds (VOCs) conducted at the Montseny Natural Park, a Mediterranean forest near Barcelona, during the summers of 2021–2023. Using Proton Transfer Reaction Mass Spectrometry (PTR-MS) and Gas Chromatography-Mass Spectrometry (GC-MS), we detected unexpected early-morning toluene peaks occurring exclusively during heatwave conditions. These peaks coincide with elevated temperatures and high vapour pressure deficits. Toluene, commonly thought to be an anthropogenically emitted volatile organic compounds, has a high potential towards tropospheric ozone formation and secondary organic aerosol formation. Therefore, we were very interested to figure out the origin of such morning peaks. We ruled out the possibility of being originated elsewhere and brought to the site due to meteorology, we also ruled out the possibility of being mainly of anthropogenic origin, which lead to the possibility that such peak is biogenic, suggesting a connection to vegetation drought stress—potentially linked to Quercus ilex emissions. Our findings provide evidence for a previously unrecognized biogenic source of toluene, with implications for atmospheric chemistry, ozone formation, and forest ecosystem resilience. We propose that early-morning toluene peaks could serve as indicators of Mediterranean forest vulnerability to climate-induced stressors.

1. Introduction

Globally, biogenic sources emit greater quantities of volatile organic compounds (VOCs) than anthropogenic sources.¹ However, both types of emissions are critical for understanding regional atmospheric chemistry. Biogenic volatile organic compounds (BVOCs), released by plants, serve diverse functions, ranging from cellular protection and defense² to chemical signaling at regional scales.³ Both BVOCs and anthropogenic VOCs (AVOCs) play a critical role in shaping interactions between the biosphere and atmosphere. These highly reactive compounds significantly influence atmospheric chemistry by serving as precursors to secondary organic aerosols (SOAs) and tropospheric ozone (O₃).

VOCs react readily with atmospheric oxidants, hydroxyl radicals (*OH) during the day, ozone throughout the day and night, and nitrate radicals (*NO₃) at night.⁴ A major consequence of these reactions is the formation of organic radicals (*RO, *RO₂) that are able to oxidise nitrogen oxide (NO) to nitrogen dioxide (NO₂). The subsequent photochemical dissociation of NO₂ yields to the formation of activated oxygen (O*), and this to the formation of tropospheric O₃ by interacting with molecular oxygen (O₂).^5–7 Exposure to high ground-level O₃ concentrations has been associated with an increase in respiratory morbidity and respiratory and cardiovascular premature mortality in humans,⁸ and with phytotoxic effects.⁹ Furthermore, O₃ plays a significant role as a contributor to global warming.¹⁰ According to the World Health Organization (WHO), both long-term exposure to O₃ levels above 60 µg m⁻³ (daily maximum 8 hours moving averages during peak 6 months period), and short-term exposure to O₃ levels above 100 µg m⁻³ (daily maximum 8 hours moving average) increase the risk of asthma attacks and hospital admissions.¹¹ Additionally, through gas-phase oxidation, VOCs can induce new particle formation via nucleation¹² and contribute to increase SOAs. In addition, VOC's oxidized products may condense onto existing particles, participate in heterogeneous reactions, be processed within clouds, or undergo further atmospheric transformation and deposition.¹³ These processes contribute to aerosol growth and cloud formation, thus affecting the Earth's radiative balance^14,15 directly by scattering and absorbing sunlight,^14,16 and indirectly, by modifying cloud formation and precipitation dynamics.¹⁷ They can also enhance primary productivity in some forest ecosystems by up to 25% through diffuse radiative fertilization.¹⁸

Among VOC's, toluene (C₆H₅CH₃) is a widely studied aromatic hydrocarbon with a high O₃ (ref. 19) and SOA forming potential under high NO_x conditions.²⁰ Toluene is traditionally considered of anthropogenic origin, with major emissions from vehicular traffic, production and use of solvents and paints, and from gasoline.²¹ Additionally, biomass burning serves as a source of toluene, particularly in the Southern Hemisphere.²² However, research suggests that toluene can also be emitted by vegetation, particularly under environmental stress conditions such as elevated temperatures and drought.^22,23 This could be particularly relevant for Mediterranean regions, where vegetation is increasingly exposed to extreme heat and water scarcity due to climate change.

The Western Mediterranean Basin is prone to elevated O₃ episodes due to low summer rainfall, intense solar radiation, high NO₂ pollution, increased BVOCs emissions, vertical recirculation of air masses driven by the region's complex orography and specific summer meteorological dynamics.^24–29 In fact, North-eastern Spain suffers from frequent O₃ extreme episodes exceeding the EU hourly threshold of 180 µg m⁻³.³⁰ Montseny Natural Park (60 km northeast of Barcelona) has been shown to have an O₃ production that increases with elevation.^31,32 This O₃ has already caused documented plant damage^31,33,34 and frequently exceeds standard regulatory thresholds.³⁵

While prior studies^22,23 have documented toluene emissions by vegetation, long-term observational studies are lacking. This work aims to elucidate the factors controlling toluene emissions from vegetation as stress indicators in the Mediterranean forest of the Montseny Natural Park during heatwave conditions of the summers, from 2021 to 2023. Especially during 2023, these ecosystems experienced a severe drought, the most intense in the last 14 years,³⁶ potentially exacerbating stress. We hypothesize that early morning increases in the toluene-to-benzene ratio during heatwaves serve as a more robust indicator of vegetation stress than toluene mixing ratios alone.

2. Methodology

2.1 Sampling location

The Montseny Natural Park, situated approximately 40 kilometres from the Mediterranean coast and about 60 kilometres northeast of the city of Barcelona (Catalonia, Spain), stands out for its diverse flora and fauna and varied topography, making it an ideal site for investigating the ecological impacts of climate variability on Mediterranean forest ecosystems. The area is characterized by a holm oak (Quercus ilex) forest, representative of Mediterranean montane ecosystems,³⁷ and is a characteristic Mediterranean rural site with a typical Mediterranean climate, including warm and dry summers, typically wet springs and autumns, mild winters, and highly irregular annual precipitation across years.³⁸ The park's ecological diversity and exposure to both local biogenic and regional anthropogenic emissions^39,40 make it well-suited for atmospheric chemistry studies examining VOC sources and transformations. The chosen sampling station (El Vilar de la Castanya, 41°46′45.63″ N, 02°21′28.92″ E, 720 m a.s.l.), referred to as MSY, lies at an elevation of 720 m above sea level (Fig. 1) and serves as a rural background monitoring site, providing valuable insights into regional background levels of atmospheric VOCs and other atmospheric pollutants.^41,42 A satellite view containing major cities and highways can be found in Fig. S.1. This site is included in the European Aerosols, Clouds, and Trace Gases Research Infrastructure Network (ACTRIS), the Global Atmosphere Watch (GAW) network of the World Meteorological Organization (WMO) and the Air Quality Monitoring Network (AQMN) of the Catalan Government.

Fig. 1 Map location and topographic profiles of Barcelona (BCN; urban background) and Montseny (MSY; regional background) measurement sites. Images credits: ©Google Earth Pro.

The MSY site has been demonstrated to be representative of the regional background and is sufficiently elevated and distant from specific urban anthropogenic emission sources to be directly affected by them.^41–48 However, the site can occasionally be affected by pollution transported from urban and industrial areas under specific meteorological conditions such as anticyclonic patterns and regional recirculation events.^{24,26,29,46,49,50} In addition, infrequent emissions from nearby vehicles may occasionally influence local measurements.

2.2 PTR-MS sampling

In this study, we conducted continuous, real-time monitoring of VOC mixing ratios at MSY using a Proton Transfer Reaction quadrupole mass spectrometer (PTR-MS, Ionicon PTR-MS Quad, Ionicon Analytik, Innsbruck, Austria), during the summer period (June to August) of 2021–2023. This technique is based on the protonation of VOC molecules by hydronium ions (H₃O⁺), which are generated in a hollow cathode ion source at high density. The VOCs in an air sample are then protonated in reactions with the hydronium ions in a drift tube (600 V, 60 °C) under the influence of an electric field (E/N = 120 Td), producing protonated VOCs and water clusters as by-products.⁵¹ The instrument was housed in a former farmhouse within Montseny Natural Park. Ambient air was drawn through a heated (40 °C) and insulated inlet tube (1/4″ inch OD) made of PFA (perfluoroalkoxy) positioned outside a window, approximately 3 meters above ground level, to ensure a representative sampling of the surrounding forest atmosphere. Blank measurements were conducted every 3 hours by passing inlet air through heated platinum pellets to correct background signals. A cylinder containing a mixture of volatile organic compounds (VOCs) was used for calibration. This custom-made standard from Apel-Riemer Environmental Inc. (Miami, USA) included formaldehyde, methanol, acetonitrile, acetaldehyde, acetone, dimethyl sulfide, isoprene, methyl vinyl ketone, methyl ethyl ketone, methyl acetate, benzene, toluene, xylenes, chlorobenzene, 3-octen-ol, p-cymene, monoterpenes and sesquiterpenes (500 ppbv). A second standard from Restek Corporation (Bellefonte, PA, USA) contained benzene, toluene, styrene, xylenes, chlorobenzene, trimethylbenzene, dichlorobenzene and trichlorobenzene (1 ppmv). These standards, included all VOCs analyzed in this study and were calibrated. The monitored compounds included methanol (m/z 33), acetonitrile (m/z 42), acetaldehyde (m/z 45), acetone (m/z 59), isoprene (m/z 69), methyl vinyl ketone and methacrolein (m/z 71), methyl ethyl ketone (m/z 73), benzene (m/z 79), toluene (m/z 93), styrene (m/z 105), xylenes (m/z 107), and monoterpenes (m/z 137).

2.3 GC-MS measurements and analysis

In addition to continuous monitoring using PTR-MS, sporadic air samples were collected and analyzed to complement and validate the real-time measurements. These samples were obtained using an automatic cartridge sampler (Sypac, France) equipped with stainless steel sorbent cartridges. Each cartridge contained 115 mg of Tenax® TA and 230 mg of SulfiCarb®, separated and secured by retaining springs and sealed with airtight caps (Markes International Inc., USA). Ambient air was pumped at a controlled flow rate of 200 mL min⁻¹ for 45 minutes. The hydrophobic nature of the adsorbents minimized water interference, ensuring sample integrity. Before sampling, the cartridges were pre-conditioned at 350 °C for 30 minutes using purified helium (100 mL min⁻¹). Blank samples were prepared by placing cartridges in the sampler without activating the pump, and toluene levels in these blanks were subtracted from the ambient air measurements. The cartridges were analyzed using an automatic thermal desorption system (TD Auto Sampler, Series 2 Ultra, Markes International Inc., Wilmington, USA) coupled with a cryotrap injector (Unity Series 2, Markes International Inc., Wilmington, USA) and a gas chromatograph-mass spectrometer (GC-MS; 7890A GC and 5975C MSD with Triple-Axis Detector, Agilent Technologies, USA).

Thermal desorption was performed in a two-stage process using a TD-Unity system (Markes International Inc.). In the first stage, VOCs were desorbed from the sample tubes at 320 °C for 30 minutes, with a flow path temperature of 200 °C, a carrier gas flow of 50 mL min⁻¹, and a split flow of 2 mL min⁻¹. Prior to desorption, the tubes underwent a pre-purge, followed by dry purging for 2 min to eliminate moisture. The minimum carrier pressure was maintained at 0.5 kPa, and the standby split was set at 25 mL min⁻¹. In the second stage, analytes were focused on a cryogenic trap, which was pre-purged for 2 min, cooled to −25 °C, and then rapidly heated at 40 °C s⁻¹ to 320 °C, where it was held for 7 min with a split flow of 5 mL min⁻¹. The desorbed compounds were then transferred to the GC column (HP 5MS, 30 m × 0.25 µm × 0.25 mm, Agilent Technologies) via a transfer line heated to 250 °C.

The GC column temperature program began with an initial hold at 35 °C for 3 minutes, followed by a ramp of 15 °C min⁻¹ to 150 °C, 50 °C min⁻¹ to 250 °C (5 minutes hold), and 30 °C min⁻¹ to 280 °C (5 minutes hold). The helium carrier gas flow rate was maintained at 1 mL min⁻¹, with a total run time of 30 minutes.

GC-MS chromatograms were analyzed using relevant databases and quantified by comparing to the same gas standard as used for the PTR-MS analysis (Apel-Riemer Environmental, Inc.).

2.4 Meteorological parameters and ancillary pollutants

Meteorological variables were recorded using a DAVIS Vantage Pro2 (Davis Instruments Corp., Hayward, California, USA) automatic weather station installed at the MSY sampling site. The weather station provided measurements of air temperature (°C), relative humidity (%), wind speed (m s⁻¹), wind direction (°), solar radiation (W m⁻²), UV index (unitless) and atmospheric pressure (bar), with a time resolution of 5 minutes. Ozone (O₃) was measured using a 48 UV analyzer MCV (Monitor for Chemical Vapors, Barcelona, Spain), which employs ultraviolet (UV) photometry to detect O₃ based on its absorption at 254 nm. Carbon monoxide (CO) concentrations were measured with a model 300EU UV photometric analyzer (Teledyne API, San Diego, CA, USA), which employs gas filter correlation infrared absorption to detect CO with high sensitivity and precision. Nitrogen oxides (NO and NO₂, jointly known as NO_x) were measured using a Thermo Scientific 42i-TL chemiluminescence analyzer, designed for trace-level detection.

In addition, equivalent black carbon (eBC) mass concentrations were monitored using the AE33 dual-spot aethalometer (Magee Scientific, Berkeley, CA, USA), which determines the attenuation of light through aerosol-loaded filter tape and provides particle light absorption coefficients at seven different wavelengths (370, 470, 520, 590, 660, 880, and 950 nm).⁵² In this study, all eBC observations were performed with the AE33 at its original 1 minute time resolution, ensuring high temporal accuracy and minimal measurement bias.

Barcelona's urban emissions can significantly influence air quality in Montseny, especially under certain meteorological conditions such as sea and mountain breeze circulation. These air masses can transport pollutants, including VOCs and other urban-derived compounds, from the city to the forested area. For this reason, in our study we focus on the boundary layer height (BLH) over Barcelona rather than Montseny. The Barcelona city BLH was obtained from the Copernicus Climate Change Service, Climate Data Store.⁵³ The specific coordinates used were: 41°23′–41°24′–N, 2°06′–2°07′–E.

2.5 Random Forest analysis

The Random Forest model was employed to investigate the relationships between environmental variables and toluene mixing ratios measured at the Montseny station. In this analysis, the dependent variable was the toluene mixing ratios, while the independent variables included meteorological factors such as temperature, relative humidity, vapor pressure deficit, wind speed, and solar radiation. One of the key variables incorporated in the Random Forest model is the Temperature–Humidity Index (THI, dimensionless index), which integrates the combined effects of air temperature (T, measured in degrees Celsius, °C) and relative humidity (RH, expressed as a percentage, %) to provide a more realistic measure of environmental stress.⁵⁴ The THI is calculated as follows:

THI = 0.8 × temp + RH × (temp − 14.4) + 46.4

To explore the relationship between predictor variables and the toluene mixing ratios we applied a Random Forest (RF) classification approach using the Random Forest package in R.⁵⁵ Prior to model development, the dataset was randomly divided into training (80%) and testing (20%) subsets using the caret package,⁵⁶ with stratification to preserve the distribution of the response variable across subsets. A fixed random seed was set to ensure reproducibility. Although Random Forest models do not rely on assumptions such as linearity or normally distributed residuals, model performance was evaluated on the independent test subset obtained from the 80–20% split.

The RF model was configured to construct 500 decision trees (n.tree = 500). The number of variables randomly sampled at each split (mtry) was set to the square root of the total number of predictor variables, following Random Forest function defaults. The model was run on the training set and evaluated on the testing set using appropriate performance metrics (e.g., accuracy, RMSE, R²), depending on the task.⁵⁵

To assess the contribution of each variable to the model's performance, we computed two measures of variable importance, (i) the permutation importance (mean decrease in accuracy), obtained by randomly permuting the values of a given predictor across observations and measuring the reduction in out-of-bag prediction accuracy, and (ii) the impurity-based importance (mean decrease in Gini), calculated as the total decrease in node impurity (Gini index) attributed to splits on that predictor, averaged across all trees in the forest. Permutation importance was computed using the default implementation in the randomForest package (type = 1), which permutes each variable across cases while keeping other predictors fixed. Model performance was evaluated using an 80–20% train–test split with stratification (see above), and by inspecting residuals and variable importance plots to verify the model's predictive reliability.

2.6 Statistical analysis

The Mann–Whitney U test, a non-parametric test for independent samples, was employed to assess differences in VOC diel median mixing ratios and other meteorological parameters and ancillary pollutants between days exhibiting early morning peaks and days without early morning peaks. This analysis was performed using the R software following the function Wilcox.test. Significant differences were determined using a threshold of p < 0.01.

3. Results

3.1 Validation of toluene measurements

One of the limitations of the PTR-MS technology is that it is mass-selective rather than compound-selective, and thus other compounds (such as fragments or isomers) may contribute to the monitored signal. The m/z 93 signal, typically assigned to protonated toluene, however other compounds and clusters can also be found at this m/z, such as propionic acid water cluster, chloroacetone, 1,2,3-propanetriol, or an ethanol cluster. More importantly m/z 93 can also receive minor contributions from the fragmentation of monoterpenes (e.g., α-pinene, 3-carene) at a typical rate of 1–2% of the monoterpene signal (Misztal et al., 2012; Ambrose et al., 2010). This same issue can be found for ρ-cymene, fenchone and thujone,^57–60 which fragmentation importance for all these has been reported for holm oaks.⁶¹ However, no explicit correction for monoterpene fragmentation was applied to the m/z 93 signal in our study. This decision is supported by parallel GC-MS validation (Fig. 2), where occasional offline GC-MS measurements were conducted at the same location and inlet as the PTR-MS system, and which independently quantified toluene and confirmed its dominance at m/z 93, demonstrating that monoterpene contributions were negligible under the conditions observed. The 1 : 1 plot of the toluene concentrations obtained using both techniques is depicted in Fig. 2. An orthogonal distance regression including uncertainties in both variables was performed to quantify the relationship between the PTR-MS toluene results and the GC-MS toluene results. The fit yielded a slope of 1.17 ± 0.21 and an intercept of −0.009 ± 0.037, indicating that the intercept is statistically indistinguishable from zero. This suggests that the two variables follow an approximately proportional relationship. The reduced chi-squared value (∼0.86) indicates that the model adequately describes the data within the assumed measurement uncertainties. Overall, the results support a strong linear correlation between the two analytical techniques. Furthermore, some GC-MS sampling coincided with early morning toluene peak days. In general, Fig. 3 shows that the time series of toluene mixing ratios measured by PTR-MS and GC-MS were generally in agreement, including the times when the early morning toluene peak was present, validating the PTR-MS measurements. These results indicate that at least for the dates for which GC-MS measurements are available, the reported mixing ratios for m/z 93 measured at MSY correspond to toluene.

Fig. 2 Relationship between PTR-MS and GC-MS. Points represent individual measurements; horizontal and vertical error bars indicate the uncertainty in of the GC-MS data (30%) and the standard deviation of PTR-MS 45 min averages (coinciding with the same 45 minutes period of the cartridge sampling), respectively. The solid red line shows the orthogonal distance regression (ODR) fit accounting for uncertainties in both variables and the black dashed line represents the 1 : 1 line. ODR results in a slope of 1.17 ± 0.21 and an intercept of −0.009 ± 0.037, with a chi-squared value of 0.86.

Fig. 3 Time series of toluene mixing ratios (ppbv) measured by PTR-MS average, in red when there was no early morning toluene peak and in purple when there was an early morning toluene peak, with standard deviations as red shade, together with the GC-MS observations of toluene in black squares with its 30% error for 2021 (top graph) and 2022 (bottom graph).

3.2 Toluene diel variability: evidence for an early morning peak

The diel cycle of toluene mixing ratios (Fig. 3) was calculated from the whole time series for the monitored period (from June 1st to August 31st, 2021–2023) (Fig. S.2.1 and S.2.2). During the 276 days observation period, 52 days had instrument malfunctions. Of the remaining days, 152 showed no toluene peak, while 72 days displayed an early morning toluene peak. Ten of these 72 peak days occurred in back-to-back periods. Furthermore, 143 days were classified as heat wave days (average daily temperature > 22 °C for >3 consecutive days), with 57 of these also exhibiting an early morning toluene peak (SI Fig. S.2.2). Toluene mixing ratios start to increase early in the morning, from approximately 07:00 local time (LT), with a slight peak during early morning and reaching their maximum between 12:00 and 15:00 LT. Afterwards, mixing ratios remain relatively stable for a few hours before showing a transient increase at 20:30 LT to rapidly decline after sunset, around 21:30 LT.

The red lines represent the 25th, 50th (median), and 75th percentiles of the hourly toluene mixing ratios. Blue triangles correspond to the average mixing ratios for each hour. A detailed inspection of the dataset revealed that some days presented very strong early morning peaks of toluene mixing ratios (from 7:00 to 9:00 LT).

3.3 Separation of peak days from non-peak days

A detailed inspection of the dataset revealed that some days presented very strong early morning peaks of toluene mixing ratios (from 7:00 to 9:00 LT). Examples of this early morning peak are shown in Fig. S.3.1 and S.3.2. Peak days were identified based on the visual presence of a distinct early morning toluene peak, characterized by a rapid increase in concentration occurring consistently before the onset of the sea breeze (typically around 09:00 LT). This approach, rather than a fixed numerical threshold, was chosen due to the inherent variability in peak magnitude influenced by other atmospheric factors. The observed early morning peaks were consistently orders of magnitude above the instrument's minimum detection limits and clearly discernible from instrumental noise, ensuring the robustness of the peak vs. non-peak classification. Limit of detection for all the VOCs are shown in Table S.1.

To better represent the typical diurnal variation and reduce the influence of extreme values, we report the median mixing ratio for each time interval. Peak and non-peak diel comparisons were generated for all available VOCs, trace gases, aerosol parameters, and meteorological variables to allow a systematic multi-species comparison (see SI Section 4). This approach was intended to assess whether the early-morning toluene peaks coincide with co-varying chemical or meteorological processes that could provide insight into their origin, rather than to characterize the full diurnal cycles of individual species. These comparisons are depicted in Fig. 4 and 5, with full boxplots provided in SI Fig. S.4.1 and S.4.2. Toluene shows a distinct early morning peak around 07:00–09:00 LT, with a maximum at 08:30 LT on peak days (Fig. 5a). The mixing ratios during early morning peak days are statistically higher (p < 0.01), with a median value of 0.23 ppbv and interquartile range (IQR) of 0.18 ppbv, compared to 0.08 (IQR 0.07) ppbv on non-peak days. A second peak is found between 13:00 and 14:00 LT with a median value of 0.26 (IQR 0.17) ppbv for peak days, compared to 0.17 (IQR 0.04) ppbv on non-peak days. A third peak, with a maximum around 20:30 LT on non-peak days, reaches a median of 0.16 (IQR 0.13) ppbv, while on peak days the median value is slightly higher at 0.21 (IQR 0.20) ppbv.

Fig. 4 Box plot showing the diel variation of toluene (m/z 93) mixing ratios (ppbv) in local time, measured by PTR-MS, during the period from June 1st to August 31st, 2021–2023. Triangles represent the mean, the dots represent the medians, the box the interquartile range and the error bars are the 75 and 25 percentile.

Fig. 5 Median diel variation of all monitored VOCs for the days that we observe an early morning toluene peak (blue) and the days when we don't observe the early morning peak in toluene (black). Units are in ppbv. (a) Toluene, (b) monoterpenes, (c) benzene, (d) methanol, (e) acetonitrile, (f) acetaldehyde, (g) acetone, (h) isoprene, (i) MACR + MVK and (j) MEK. Shaded areas represent the 75th and 25th percentile. Solid black line represents the LOD for each compound.

Fig. 5b presents the diurnal cycle of monoterpenes, revealing a distinct early morning peak on peak days between 07:30 and 10:00 LT, closely resembling the pattern observed for toluene. This peak reaches a maximum around 08:30–09:00 LT, with median values of 0.96 (IQR 1.18) ppbv, compared to a significantly lower median of 0.22 (IQR 0.23) ppbv on non-peak days (p < 0.001). The early morning toluene peaks often mirrored the diurnal pattern of monoterpenes. However, it is important to note that early morning monoterpene peaks were frequently observed even on days without an early morning toluene peak. This further supports that the PTR-MS m/z 93 signal, while potentially having a minor monoterpene contribution, is predominantly driven by toluene and not significantly compromised by monoterpene interference.

Benzene (Fig. 5c) showed elevated mixing ratios on toluene peak days, with mixing ratios gradually increasing in the morning (07:00–10:30 LT) and reaching a maximum between 14:00–14:30 LT. On peak days, the median mixing ratio of benzene was 0.05 (IQR 0.04) ppbv, while on non-peak days it was slightly lower at 0.05 (IQR 0.03) ppbv (p < 0.01).

Methanol (Fig. 5d) exhibits a pronounced diurnal pattern. During the night and early morning hours (00:00–07:00 LT), mixing ratios remain low, particularly on non-peak days. From 07:00 LT onward, methanol levels begin to rise, with a distinct morning peak occurring between 08:30 and 09:00 LT. During this interval, peak days display significantly higher mixing ratios of 3.99 (IQR 3.42) ppbv compared to non-peak days of methanol 3.19 (IQR 1.92) ppbv (p < 0.001). Methanol levels continue to increase after the morning hours, reaching a broader midday maximum around 13:30 LT.

Acetonitrile (Fig. 5e) shows a clear diurnal variation, with mixing ratios beginning to rise after 06:00 LT on peak days. A morning peak occurs between 06:30 and 10:30 LT, reaching its maximum around 08:00 LT. During this period, peak days exhibit significantly higher mixing ratios of acetonitrile, 0.14 (IQR 0.05) ppbv, compared to non-peak days, 0.12 (IQR 0.04) ppbv (p < 0.001), suggesting the presence of episodic sources such as biomass burning events, or other irregular activities that can temporarily elevate background levels of acetonitrile. A broader secondary maximum is also observed around midday at 13:00 LT. Notably, a detailed inspection (Fig. S.4.1k and l) shows that the observed mean in acetonitrile morning acetonitrile peaks were predominantly driven by two specific episodes, one in 2021 and another in 2022, suggesting that biomass burning or anthropogenic events may have contributed disproportionately to the average peak-day signal.

Acetaldehyde (Fig. 5f) and acetone (Fig. 5g) follow similar diel profiles, with early morning increases, midday peaks (around 13:30 LT for both compounds), and afternoon declines. Acetone reaches its maximum peak at 13:30 LT, with values of 5.5 (IQR 3.56) ppbv on peak days and 4.22 (IQR 3.58) ppbv on non-peak days (p < 0.01). Acetaldehyde also reaches its peak around the same time, with values of 2.21 (IQR 0.65) ppbv on peak days, higher compared to 1.71 (IQR 0.90) ppbv on non-peak days (p < 0.001).

Isoprene (Fig. 5h) shows a distinct morning increase on both peak and non-peak days, starting around 07:30 LT and continuing until 10:00 LT, with a maximum observed between 13:00–13:30 LT. The median mixing ratio on peak days is 0.61 (IQR 0.60) ppbv, while on non-peak days it is significantly lower at 0.37 (IQR 0.29) ppbv (p < 0.01).

Methacrolein and methyl vinyl ketone (MACR + MVK, Fig. 5i) and Methyl Ethyl Ketone (MEK, Fig. 5j) exhibit similar morning increases and peaks during midday. MACR + MVK reaches its maximum peak at 13:00 LT, with a median value of 0.82 (IQR 0.96) ppbv on peak days, higher compared to 0.50 (IQR 0.39) ppbv on non-peak days (p < 0.001). MEK presents an early increase between 07:30 and 11:30 LT, followed by a broader midday maximum, peaking at 14:30 LT. On peak days, median values reach 0.65 (IQR 0.05) ppbv, higher compared to 0.47 (IQR 0.35) ppbv on non-peak days (p < 0.01).

3.4 Comparison of meteorological conditions and ancillary pollutants on days with and without morning toluene peaks

The diurnal variation of temperature (Fig. 6a), vapor pressure deficit (VPD) (Fig. 6b), wind speed (Fig. 6c), wind direction (Fig. 6d), NO (Fig. 6e), NO₂ (Fig. 6f), NO_x (Fig. 6g), black carbon (eBC) (Fig. 6h), CO (Fig. 6i), ozone (Fig. 6j), radiation (Fig. 6k) and BLH (Fig. 6l) during days with and without early morning toluene peaks are depicted in Fig. 6. The temperature box plots showed a typical diurnal cycle, with values gradually increasing after sunrise, peaking in the early afternoon (13:30 LT), and decreasing during the evening and night. On peak days, the temperature appeared higher (30 IQR 5.18 °C, at 13:30 LT) than on non-peak days (25 IQR 5.11 °C, at 13:30 LT) (Fig. 6a) (p < 0.001). The atmospheric vapor pressure deficit (VPD), calculated from the combination of relative humidity and temperature, followed a similar diurnal pattern, with values rising throughout the morning and reaching their peak in the afternoon before gradually declining. The peak-day VPD was generally higher compared to non-peak days (p < 0.001), especially between midday and the late afternoon maximum (peak days 2526 (IQR 1533) Pa for peak days and 1649 (IQR 974) Pa for non-peak days) (Fig. 6b). The wind direction plots indicated relatively stable north conditions during the night and early morning. The differences shown for the peak days are simply due to a change from north to north-northeast. Then for both peak and non-peak conditions wind direction come from south-southwest while sea and mountain breezes are developed. Also, there was a noticeable shift in wind speed after 09:00, with peak and non-peak days exhibiting similar trends, although slightly higher wind speeds could be found on peak days around 14:00–15:00 LT (p < 0.01). During the night and early morning hours (00:00–06:00 LT), BC concentrations remained relatively stable around 0.3–0.4 µg m⁻³. After 06:00 LT, a gradual increase was observed reaching a maximum around 13:00 LT, followed by a decline towards the evening. BC mean levels appeared slightly higher on peak days with 0.57 (IQR 0.31) µg m⁻³ compared to the 0.45 (IQR 0.25) µg m⁻³ for non-peak days (Fig. 6h) (p < 0.01).

Fig. 6 Median diel variation of (a) air temperature (in °C), VPD in Pa (b), wind speed in m s⁻¹ (c), wind direction in ° (d), NO in µg m⁻³ (e), NO₂ in µg m⁻³ (f), NO_x in µg m⁻³ (g), black carbon (BC) in µg m⁻³ (h), CO in µg m⁻³ (i), ozone in µg m⁻³ (j), radiation in W m⁻² (k) and BLH in m (l), measured at Montseny during days with (blue) and without toluene morning peaks (black) (note that BLH is from Barcelona city). Shaded areas represent the 75th and 25th percentile. Note that for wind direction only averages and standard deviation is expressed.

The ozone diel cycle showed lower concentrations during the night and early morning, followed by a steady increase throughout the day, reaching peak levels at around 15:00 LT. The concentrations were higher for peak days, 112.5 (IQR 32.4) µg m⁻³, than for non-peak days, 109.9 (IQR 34.6) µg m⁻³ (Fig. 6j) (p < 0.01). In order to assess the ozone exposure diel variations of 8 h moving averages of O₃ for peak and non-peak days are depicted in Fig. S.6. The maximum 8 h moving averages (MDA8) are higher on peak days (Fig. S.6). NO, NO₂ and NO_x concentrations were higher during the evening hours, particularly between 12:00–18:00 and 18:00–23:00. On non-peak days, NO_x showed slightly elevated concentrations 3.91 (IQR 3.16) µg m⁻³ compared to peak days (Fig. 6g) 3.18 (IQR 2.90) µg m⁻³ (p < 0.01), particularly during the late afternoon and evening hours. Similarly, NO₂ concentrations were slightly higher during the same period on non-peak days 3 (IQR 2.17) µg m⁻³, compared to the peak days (2.93 (IQR 2.09) µg m⁻³) (Fig. 6f) (no statistically significant difference). No apparent difference for NO was found between peak and non-peak day (Fig. 6e) (no statistically significant difference). During days with a morning toluene peak, solar radiation values (Fig. 6k) are generally higher throughout the day compared to days without a morning peak. The difference is most pronounced during midday (12:00–15:00 LT), radiation reaches approximately 900–950 W m⁻² on peak days, while on non-peak days it remains lower, around 800–850 W m⁻². This pattern suggests that higher incoming solar radiation in the morning may contribute to enhanced photochemical activity or increased plant emissions, potentially triggering the observed toluene morning peaks, nevertheless this difference in the medians is not statistically significant.

The BLH in Barcelona appeared not to have a statistically significant difference in the overall diurnal evolution of BLH between days with a morning toluene peak and days without it (Fig. 6l). In both cases, BLH started low during the night, began to rise around 06:00, reached a maximum during midday, between 12:00 and 15:00 and it gradually decreased toward the evening. The median height for days where the toluene early morning peak was present was 661 m whereas the days with no toluene early morning peak the BLH was 633 m (p > 0.005). In any case, at 8:30 LT on peak days, the BLH from Barcelona seemed lower (343 m) than the MSY station elevation (720 m a.s.l.).

The toluene to benzene ratio (T/B) was calculated (Fig. 7). During peak days, a pronounced increase in the T/B ratio was observed in the early morning hours (7:00–10:00 LT), with median values reaching up to 6.5. After this morning's enhancement, the T/B ratio gradually returns to values similar to those of non-peak days. On non-peak days, the T/B ratio remained relatively constant between 3.3 and 4.

Fig. 7 Median diel variation of the toluene-to-benzene (T/B) ratio during peak (blue) and non-peak (black) days.

3.5 Drivers of morning toluene peaks, results from Random Forest

The Temperature–Humidity Index (THI) and the variables of cumulative temperature and relative humidity sums (temp_sum_6h, RH_sum_62h, temp_sum_48h) were identified as the most influential predictors, showing the highest values for both %IncMSE and IncNodePurity (Fig. 8). According to the variable importance analysis from the Random Forest model. This indicates that the plant response integrates thermal and drought stress over time rather than reacting only to instantaneous conditions. This pattern supports the idea that the magnitude of the toluene peak is related to the duration and persistence of the preceding heatwave.

Fig. 8 Variable importance in Random Forest model: %IncMSE (in salmon) and IncNodePurity (in blue). Summer dataset, the model explains 81.9% of the variance of toluene peaks.

While environmental stressors like ozone exposure are known to influence biogenic emissions, incorporating metrics such as previous-day peak O₃ or MDA8 O₃ into our hourly Random Forest model presented methodological challenges. Directly integrating daily aggregated O₃ data into an hourly model would have introduced a mismatch in temporal resolution and potential temporal leakage. While future work could explore lagged O₃ effects with a re-architected modeling framework, it was deemed outside the scope of the current hourly prediction model. The insights from same-day MDA8 O₃ (SI Fig. S.6) primarily served as descriptive context rather than a direct predictor in our RF analysis.

4. Discussion

Toluene is a widely studied aromatic hydrocarbon, traditionally regarded as an anthropogenic compound, with its atmospheric presence largely attributed to human activities. Consequently, toluene is often used as an indicator of anthropogenic influence and industrial activity, providing valuable insight into air quality. While these human-related emissions are well documented, increasing evidence supports the existence of biogenic toluene emissions.^22,23 Indeed, toluene emissions have already been reported for several plant species, including Helianthus annuus (sunflower),²² Quercus ilex (holm oak),⁶¹ Halimium halimifolium (spotted yellow sun rose),⁶² Thymus vulgaris (thyme),⁶³ Triticum aestivum L. (common wheat),⁶⁴ and toluene increases have been observed during summer months in rural ecosystems partly attributed to increased biogenic emissions.⁶⁵ Therefore, the possibility of a biogenic origin of the toluene peaks observed during early mornings on the days of heat waves is plausible. It is important to note, however, that the measurements in this study represent atmospheric mixing ratios, which are influenced not only by emission of VOCs (both anthropogenic and biogenic) and deposition, but also by atmospheric mixing and chemical transformations driven by meteorological conditions and oxidant levels are certainly affecting the observed dynamics.⁶⁶

4.1 Local vs. regional origin of toluene peaks

Air circulation at MSY is strongly influenced by the valley's topography. During active mountain and sea breeze hours (starting at approximately 9:00 LT), VOCs, aerosols and trace gases are mainly transported from the urban and industrial areas of Barcelona.^38,39 This is evident in the diel cycles of toluene, BC, NO_x, NO, CO, and O₃—compounds typically associated with pollution—as well as observable in other anthropogenic VOCs such as benzene, acetonitrile, acetaldehyde, acetone, MVK + MACR, MEK, styrene and xylenes. On peak days, temperature, VPD, and wind speed were higher than on non-peaks days, reflecting heatwave conditions. Higher temperatures enhance convection, leading to stronger winds and boundary layer mixing.

The dynamics of the planetary boundary layer, which varies in height throughout the day, play a crucial role in the dilution and accumulation of emitted compounds.^67,68 During night and early morning, a shallow boundary layer can trap emissions near the surface,⁶⁹ resulting in higher mixing ratios that later decrease as the boundary layer grows. During this study, analysis of the BLH over Barcelona urban area revealed that, at the time of the early morning toluene peaks, the BLH was below the altitude of MSY. This, together the sea breeze, which typically initiates around 09:00 LT, plays a critical role in pollutant dispersion and vertical mixing at the site,²⁷ indicates that polluted air masses from Barcelona were likely trapped below and not transported upward to MSY. Thus, the observed early morning toluene peaks are most likely of local origin.

Further evidence comes from the toluene-to-benzene (T/B) ratio, a common indicator of air-mass photochemical age. Benzene, being relatively stable, has a longer atmospheric lifetime than toluene, which reacts more rapidly with hydroxyl radicals.^4,70 Freshly emitted air masses exhibit high T/B ratios, whereas aged air shows lower ratios due to toluene oxidation.⁷¹

Previous studies reported a T/B ratio of 3.58 at MSY,⁴⁸ consistent with our observed range (3.3–4). However, during early morning peaks, the ratio increased to a median of 6.52, suggesting either local toluene emissions or temporarily reduced OH concentrations that slow toluene oxidation.

4.2 Impact of atmospheric chemistry in toluene dynamics

In addition to the early morning peak, two other toluene peaks were observed in the diel cycle: one around midday, associated with transport of polluted air from Barcelona, and another around 20:00–21:00 LT, visible also in monoterpenes, benzene, isoprene, methanol, and MVK + MACR. This late-evening peak coincided with increases in NO, NO₂, and consequently NO_x.

Two processes may explain this pattern:¹ a reduction in photooxidation as sunlight declines, and² a reduction in atmospheric dilution due to the evening decrease in BLH. However, this peak was not observed for acetonitrile, MEK, or acetone, suggesting that reduced dilution alone does not fully explain it.

Since OH is mainly produced by ozone and water photolysis under light, lower evening radiation reduces OH and O₃ availability, thereby slowing VOC oxidation. This can lead to temporary accumulation of reactive VOCs such as toluene, monoterpenes, isoprene, and MVK + MACR. The higher T/B ratio observed in the late evening further supports a reduction in OH concentrations. Additionally, it is important to understand the non-linear photochemical behaviour of toluene, where under low NO_x, toluene oxidation by OH scavenging is further reduced by compounds such as isoprene.⁷²

The early morning peak, in contrast, likely reflects a combination of reduced photochemistry and strong local emissions. Both toluene and monoterpenes showed synchronized early morning peaks, which rapidly declined once photochemical activity increased. Afterwards the transport of polluted air masses by sea and mountain breeze adds to the mixing ratios of the not so reactive monoterpenes and toluene, together with a possible continued local emission during the day. Although monoterpenes are highly reactive with ozone, toluene's lower ozone reactivity⁷³ explains why its decline was less pronounced during photochemically active hours. Additionally, night-time values remain low, flat, and show no systematic accumulation, even during shallow nocturnal boundary-layer conditions. This lack of nocturnal build-up is an indication that emissions are minimal at night and that the observed peaks are not caused by storage and nocturnal accumulation. Instead, the sharp early morning increases suggest the involvement of emission processes that are active during daylight hours, likely linked to light, temperature, and physiological responses of the vegetation.

4.3 Anthropogenic vs. biogenic origin of toluene peaks

Having established a local origin, we examined whether the early morning toluene peaks were anthropogenic or biogenic. Some anthropogenic tracers—such as acetonitrile, benzene, and BC—also showed early morning increases for their averages only on certain days, (see Fig. S.5) on the 26th July 2021 and the 18th July 2022 (for BC only on 18th July 2022 given the lack of data for 2021), thus driving the averages higher. This suggests that episodic events like biomass burning or localized anthropogenic activities may have temporarily elevated background levels of these compounds on those specific days. Since this effect was only seen during two days, the rest of the days with toluene early morning peak seemed to not be related to pollution events.

Conversely, similar early morning peaks were observed in monoterpenes, known biogenic compounds previously reported to be emitted by holm oaks in MSY.^{38,39,48,61,74} Although possible mass spectral interference between monoterpenes and toluene at m/z 93(ref. 57, 58 and 61) was considered, GC-MS measurements confirmed comparable magnitudes in both techniques, minimizing this concern.

These peaks were most pronounced during heatwaves, when plants experience thermal and water stress that enhance VOC emissions. Interestingly, monoterpene peaks appeared at the onset of heatwaves, whereas toluene peaks emerged only after several days of sustained stress. Their timing coincided with sunrise, when photosynthesis and VOC biosynthesis begin. Methanol, another biogenic VOC linked to leaf development and cell wall expansion,^75–77 also exhibited early morning peaks, reinforcing the interpretation of a biogenic origin for toluene emissions.

Additionally, aromatic compounds emitted by plants are generally known to be synthesized within the shikimate pathway.⁷⁸ However, the divergence of toluene dynamics from other aromatic compounds suggests the shikimate pathway is unlikely the pathway for toluene biosynthesis.²² Instead, evidence suggests that toluene is likely linked to photosynthetic activity and can be produced from de novo biosynthesis and emission is not mainly driven by storage pool release.^22,23,62 It has been hypothesized and further supported that toluene is emitted in the MEP pathway, similar to monoterpenes, given the high correlation found among these two VOC types in plant laboratory experiments.^22,23,62 Thus, the similarity between monoterpenes and toluene could be derived from the fact that they are biosynthesized in the same pathway.

4.4 Possible mechanisms for the occurrence of biogenic toluene emissions at the Montseny Natural Park

Given the possible biogenic origin of the toluene early morning peaks, we investigated why this peak was only apparent during ongoing heat waves. Higher temperatures not only enhance biogenic emissions due to faster metabolic activity due to increased temperature, but also it can exert plant stress, which results in higher VOC emission as plants use VOCs to alleviate biotic and abiotic stress.² Furthermore, higher temperature can also accelerate photooxidation processes, leading to faster chemical transformation of VOCs and its degradation products in the atmosphere.⁷⁹ Furthermore, the MDA8 values are higher during peak days (where plants are exposed to higher ozone levels), thereby possibly inflicting an additional phytotoxic stress. Several studies suggest that toluene vegetation emissions are particularly enhanced under high-temperature conditions, supporting the hypothesis of toluene as a stress-induced compound.^22,23,64 Toluene emissions have been observed to increase under stress conditions in sunflowers and pines (Pinus spp.),²² as well as in Populus balsamifera in response to heat stress and to combined spider mite and subsequent heat stress.²³ Heiden et al.,²² 1999 proposed the existence of internal toluene storage pools within plants, specifically identified sunflower (Helianthus annuus) and Scots pine (Pinus sylvestris) as examples. Therefore, it is possible that the early morning toluene peaks originate from both these storage pools and de novo biosynthesis. Aromatic VOCs can also originate from soil and litter microbial activity, and several studies have shown that microbes are capable of producing toluene and related compounds under certain environmental conditions.^80,81 Therefore, a microbial contribution at our site cannot be fully excluded. However, the night-time concentrations of toluene remained consistently low, with no evidence of nocturnal accumulation despite shallow boundary layer conditions (see Section 3.2). This behaviour contrasts sharply with the distinct, light- and temperature-driven early morning peaks. While microbial activity may contribute to the background aromatic levels, the observed diel pattern strongly suggests that microbial emissions are unlikely to be the dominant driver of the early morning toluene peaks.

The Random Forest model in this study explained 81.9% of the variance of toluene, indicating a strong predictive capacity for identifying the key drivers during the summer periods. Among the most influential predictors we found THI and the variables of cumulative temperature and relative humidity sums (temp_sum_6h, RH_sum_62h, temp_sum_48h). These variables are associated with prolonged heat stress and plant water availability,²³ reinforcing the hypothesis that biogenic toluene emissions increase under heat wave stress conditions. Cumulative temperature metrics likely capture extended heat exposure that can stimulate plant metabolic responses, including enhanced VOC release.⁸² Similarly, sustained low relative humidity levels reflect increased evaporative demand, affecting plant transpiration and stomatal conductance, both of which regulate VOC emissions. However, it is important to underline that this mechanism does not apply uniformly across all VOCs. Radiation also plays a dual role by driving both photosynthesis and photochemical activity, potentially influencing both the production and degradation of toluene in the atmosphere. Although wind-related variables such as wind speed and direction were ranked lower, they still contribute to the dispersion and transport of VOCs. Overall, the Random Forest model proved to be an effective approach for non-linear relationships among environmental stressors and toluene variability, further supporting the hypothesis that forests emit toluene when stressed from high temperatures and drought.

4.5 Limitations

This study aims to elucidate the occurrence of toluene early morning peaks by first allocating a local origin, then attributing a potential biogenic source and hypothesizing a possible mechanism for such emission regarding plant thermal and water stress. Nevertheless, it is essential to understand this study monitors only ambient mixing ratios. This means that even if evidence points out to the potential biogenic toluene emission by the Montseny natural park, anthropogenic contributions cannot be excluded together with the decreased effect of dilution and reduced photochemistry regarding OH and ozone oxidation. In fact, most likely the biogenic toluene emissions do not only happen during the early morning but rather along the day, and are combined with the transported toluene from elsewhere, as well as with potential anthropogenic local sources, all subjected to OH oxidation and the dilution effect during the day. In addition, while various other benzenoid compounds are known to be biogenically emitted by vegetation, our PTR-MS measurements focused on a specific set of m/z ratios, including benzene (m/z 79), C9 aromatics (m/z 107), and styrene (m/z 105), in addition to toluene (m/z 93). Unfortunately, the mixing ratios of C9 aromatics and styrene were consistently below their respective detection limits (SI Table S.1), precluding any analysis of their diurnal patterns. For benzene, while the mean diurnal cycle in certain datasets showed minor early morning increases, these were not consistently reflected in the median diurnal cycles, and were frequently correlated with other anthropogenic tracers like acetonitrile and black carbon (SI Fig. S.4.1 and S.4.2), indicating likely anthropogenic origins rather than a distinct biogenic pattern akin to toluene. Thus, for further studies, measuring a higher suit of aromatics compounds may help towards the understanding of plant benzenoid emissions. Furthermore, more analysis is needed to corroborate further the toluene measurements with GC-MS to totally exclude any monoterpene interference as only sporadic GC-MS sampling was available for this study. Additionally branch and leave level emission measurements are need to understand how other biogenic sources such as soil or microbes could be contributing to the signal as well as to set up experiments with control, heat, drought and heat and drought treatments to test our finding. Specially, because the Random Forest model identifies important predictors but does not establish causality, therefore potential confounding factors might influence the relationship between environmental variables and toluene emissions. Within this context, if we can further prove the hypothesis stated here, of the identified role of toluene as a stress indicator, it could not only underline the ecological health of Mediterranean forests, but also highlight its potential implications for air quality regulation during increasingly frequent heat events.

5. Conclusions

This study investigates the occurrence of early-morning toluene peaks observed in the Montseny Mediterranean forest during the summer months (June, July, and August) from 2021 to 2023. The onset of morning sunlight activates plant physiological processes, leading to a sharp increase of toluene mixing ratios. This phenomenon, only visible during ongoing heatwaves, occurs when high ambient temperatures and low relative humidity intensify plant stress responses. In the Montseny region, the dominant vegetation is Quercus ilex (holm oak), a species well-adapted to the Mediterranean climatic environment. In this study, we point to the possibility that such early morning toluene mixing ratios peaks may be derived from Quercus ilex emissions under drought stress, given the local and biogenic origin of such peak. Our results also reveal that early-morning peaks in toluene are closely tied to high temperatures and low humidity, which could indicate a biogenic response due to plant stress. Later in the day, the toluene mixing ratios are mostly dominated by anthropogenic emissions being transported to the site.

The complex biochemistry and atmospheric chemistry underlying toluene dynamics highlight the need for further investigation to fully understand its sources and transformations. Understanding the exact pathways and conditions leading to toluene biogenic production will provide deeper insights into its ecological and atmospheric role, especially in Mediterranean environments like Montseny Natural Park, where vegetation stress is an increasing concern during the summer season, and such conditions are expected to increase in the near future. Further research on biogenic volatile organic compounds (BVOCs) could enhance our understanding of ecosystem dynamics and their responses to environmental conditions. This study opens new questions: could early-morning toluene emission peaks in Mediterranean forests be an indicator of forest drought stress? And, under increasing heatwaves and drought, could these peaks be silent contributors to atmospheric reactivity, giving the high SOA and ozone formation potential of toluene?⁴⁸ Given its dual anthropogenic and biogenic sources, understanding toluene's atmospheric behaviour is essential for assessing its role in air pollution, regional climate effects, and potential feedback mechanisms between vegetation and atmospheric chemistry. In any case, management practices should incorporate monitoring toluene levels in response to climate change, aiding in the development of conservation strategies for vulnerable Mediterranean ecosystems. Such efforts could inform policy-making and offer a meaningful indicator of forest health in the face of increasing climatic stress.

Author contributions

Conceptualization: AMYS, JP. Investigation: AMYS, ML, IF, RS, AA, XQ, JP, JM, JL, NP. Data curation: AMYS, NP, ML, IF, JL, MS, MP, RS. Funding acquisition and administration: AMYS, AA, XQ, JP, MPE, NM. Methodology: AMYS, ML, IF, RS, MFM, JL, JP, NP, MS, MP, RS, JM. Visualization: AMYS, ML, RS, IF. Formal analysis: AMYS, ML, JP, IF, RS, JM. Resources: AMYS, AA, XQ, JP, MPE, IJ, NM, IF. Supervision: AMYS, JP, AA, XQ, NM. Writing original draft: AMYS, ML, JP, IF, JL, RS. Writing review & editing: all authors.

Conflicts of interest

There are no conflicts to declare.

Data availability

VOC data is available at EBAS (https://ebas.nilu.no/). The rest of the data is available at https://doi.org/10.5281/zenodo.17348547.

Supplementary information (SI) is available. See DOI: https://doi.org/10.1039/d5ea00168d

Acknowledgements

This research was also supported the Agencia Estatal de Investigación from the Spanish Ministry of Science and Innovation, FEDER funds under the projects CAIAC (PID2019-108990RB-I00), the Generalitat de Catalunya (AGAUR 2021 SGR 00447), and Departament de Territori, Habitatge i Transició Ecològica (DTER) of the Generalitat de Catalunya. AMYS acknowledges La Caixa Foundation Junior Leader retaining fellowship and Grant CEX2018-000794 S funded by MCIN/AEI/10.13039/501100011033, the support from the Consolidación Investigadora project (CNS2022-135757) and her Ramon y Cajal grant (RYC2021-032519-I). She also acknowledges the AIRPHONEMA project: PID2022-142160OB-I00 funded by MCIN. JP was supported by the Spanish Government grant PID2022-140808NB-I00 funded by MCIN, AEI/10.13039/501100011033 European Union Next Generation EU/PRTR. Authors also acknowledge the data from the Air Quality Monitoring Network (AQMN) of the Catalan Government. M. F.-M. was supported by the European Research Council project ERC-StG-2022-101076740 STOIKOS, a Junior Leader fellowship from “la Caixa” Foundation (ID 100010434), codeLCF/BQ/PI21/11830010, and a Ramón y Cajal fellowship (RYC2021-031511-I) funded by the Spanish Ministry of Science and Innovation, the NextGenerationEU program of the European Union, the Spanish plan of recovery, transformation and resilience, and the Spanish Research Agency. RS acknowledges a Ramón y Cajal grant (RYC2020-029216-I) funded by MICIU/AEI/10.13039/501100011033 and by “ESF Investing in your future,” and project PID2021-122892NA-I00 funded by MICIU/AEI and by “ERDF A way of making Europe”.

References

1. J. M. Kelly, R. M. Doherty, F. M. O'Connor and G. W. Mann, The impact of biogenic, anthropogenic, and biomass burning volatile organic compound emissions on regional and seasonal variations in secondary organic aerosol, Atmos. Chem. Phys., 2018, 18(10), 7393–7422.
2. J. Laothawornkitkul, J. E. Taylor, N. D. Paul and C. N. Hewitt, Biogenic volatile organic compounds in the Earth system, New Phytol., 2009, 183(1), 27–51.
3. A. M. Yáñez-Serrano, E. Bourtsoukidis, E. G. Alves, M. Bauwens, T. Stavrakou and J. Llusià, et al., Amazonian biogenic volatile organic compounds under global change, Global Change Biol., 2020, 26(9), 4722–4751.
4. R. Atkinson, Atmospheric chemistry of VOCs and NO_x, Atmos. Environ., 2000, 34(12–14), 2063–2101.
5. U. Epa, Air Quality Criteria for Ozone and Related Photochemical Oxidants Volume II of III, Environ. Prot., 2006, I, 820.
6. EPA, Air Quality Criteria for Ozone and Related Photochemical Oxidants, Mol. Cell. Biochem., 2006, 9(3), 131–133.
7. A. Q. Criteria, P. Matter and V. Ii, Air Quality Criteria for Particulate Matter Volume II of II, Environ. Prot., 2004, II, 1148.
8. J. J. Zhang, Y. Wei and Z. Fang, Ozone pollution: a major health hazard worldwide, Front. Immunol., 2019, 10, 1–10.
9. E. Agathokleous, Z. Feng, E. Oksanen, P. Sicard, Q. Wang and C. J. Saitanis, et al., Ozone affects plant, insect, and soil microbial communities: a threat to terrestrial ecosystems and biodiversity, Sci. Adv., 2020, 6(33), 1–17.
10. G. Jia, E. Shevliakova, P. Artaxo, N. D. Noblet-Ducoudré, R. Houghton, J. House, et al., Land–climate interactions, in Climate change and land: an IPCC special report on climate change, desertification, land degradation, sustainable land management, food security, and greenhouse gas fluxes in terrestrial ecosystems, 2019.
11. J. Miles, WHO global air quality guidelines. Part matter (PM25 PM10), ozone, nitrogen dioxide, sulfur dioxide carbon monoxide, 2021, pp. 1–360.
12. J. Kirkby, J. Duplissy, K. Sengupta, C. Frege, H. Gordon and C. Williamson, et al., Ion-induced nucleation of pure biogenic particles, Nature, 2016, 533(7604), 521–526.
13. M. Hallquist, J. C. Wenger, U. Baltensperger, Y. Rudich, D. Simpson and M. Claeys, et al., The formation, properties and impact of secondary organic aerosol: current and emerging issues, Atmos. Chem. Phys., 2009, 9, 5155–5236.
14. M. Kulmala, J. Kontkanen, H. Junninen, K. Lehtipalo, H. E. Manninen and T. Nieminen, et al., Direct observations of atmospheric aerosol nucleation, Science, 2013, 339(6122), 943–946.
15. U. Pöschl, S. T. Martin, B. Sinha, Q. Chen, S. S. Gunthe and J. A. Huffman, et al., Rainforest Aerosols as Biogenic Nuclei of Clouds and Precipitation in the Amazon, Science, 2010, 329(5998), 1513–1516.
16. O. Boucher, D. Randall, P. Artaxo, C. Bretherton, G. Feingold, P. Forster, et al., Clouds and Aerosols, in Climate Change 2013: The Physical Science Basis Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC), ed. T. F. Stocker, D. Qin, G. K. Plattner, M. Tignor, S. K. Allen, J. Boschung, et al., Cambridge University Press, Cambridge, UK and New York, USA, 2013, pp. 571–657.
17. J. Haywood and O. Boucher, Estimates of the direct and indirect radiative forcing due to tropospheric aerosols: a review, Rev. Geophys., 2000, 38(4), 513–543.
18. A. Rap, C. E. Scott, C. L. Reddington, L. Mercado, R. J. Ellis and S. Garraway, et al., Enhanced global primary production by biogenic aerosol via diffuse radiation fertilization, Nat. Geosci., 2018, 11(9), 640–644.
19. M. A. Venecek, W. P. L. Carter and M. J. Kleeman, Updating the SAPRC Maximum Incremental Reactivity (MIR) scale for the United States from 1988 to 2010, J. Air Waste Manage. Assoc., 2018, 68(12), 1301–1316.
20. R. G. Derwent, M. E. Jenkin, S. R. Utembe, D. E. Shallcross, T. P. Murrells and N. R. Passant, Secondary organic aerosol formation from a large number of reactive man-made organic compounds, Sci. Total Environ., 2010, 408(16), 3374–3381.
21. DEFRA, Non-methane Volatile Organic Compounds in the UK, 2020, pp. 1–92.
22. A. C. Heiden, K. Kobel, M. Komenda, R. Koppmann, M. Shao and J. Wildt, Toluene emissions from plants, Geophys. Res. Lett., 1999, 26(9), 1283–1286.
23. P. K. Misztal, C. N. Hewitt, J. Wildt, J. D. Blande, A. S. D. Eller and S. Fares, et al., Atmospheric benzenoid emissions from plants rival those from fossil fuels, Sci. Rep., 2015, 5, 12064.
24. M. M. Millán, M. José Sanz, R. Salvador and E. Mantilla, Atmospheric dynamics and ozone cycles related to nitrogen deposition in the western Mediterranean, Environ. Pollut., 2002, 118(2), 167–186.
25. M. M. Millán, Extreme hydrometeorological events and climate change predictions in Europe, J. Hydrol., 2014, 518, 206–224.
26. M. M. Millán, R. Salvador, E. Mantilla and G. Kallos, Photooxidant dynamics in the Mediterranean basin in summer: results from European research projects, J. Geophys. Res. Atmos., 1997, 102(7), 8811–8823.
27. X. Querol, G. Gangoiti, E. Mantilla, A. Alastuey, M. C. Minguillón and F. Amato, et al., Phenomenology of high-ozone episodes in NE Spain, Atmos. Chem. Phys., 2017, 17(4), 2817–2838.
28. C. Pérez, M. Sicard, O. Jorba, A. Comerón and J. M. Baldasano, Summertime re-circulations of air pollutants over the north-eastern Iberian coast observed from systematic EARLINET lidar measurements in Barcelona, Atmos. Environ., 2004, 38(24), 3983–4000.
29. G. Gangoiti, M. M. Millán, R. Salvador and E. Mantilla, Long-range transport and re-circulation of pollutants in the western Mediterranean during the project Regional Cycles of Air Pollution in the West-Central Mediterranean Area, Atmos. Environ., 2001, 35(36), 6267–6276.
30. J. Massagué, E. Torre-Pascual, C. Carnerero, M. Escudero, A. Alastuey and M. Pandolfi, et al., Extreme ozone episodes in a major Mediterranean urban area, Atmos. Chem. Phys., 2024, 24(8), 4827–4850.
31. M. Díaz-de-Quijano, J. Peñuelas and À. Ribas, Increasing interannual and altitudinal ozone mixing ratios in the Catalan Pyrenees, Atmos. Environ., 2009, 43(38), 6049–6057.
32. À. Ribas and J. Peñuelas, Surface ozone mixing ratio increase with altitude in a transect in the Catalan Pyrenees, Atmos. Environ., 2006, 40(38), 7308–7315.
33. B. S. Gimeno, J. Peñuelas, J. L. Porcuna and R. A. Reinert, Biomonitoring ozone phytotoxicity in eastern Spain, Water, Air, Soil Pollut., 1995, 85(3), 1521–1526.
34. M. J. Sanz, V. Calatayud and E. Calvo, Spatial pattern of ozone injury in Aleppo pine related to air pollution dynamics in a coastal-mountain region of eastern Spain, Environ. Pollut., 2000, 108(2), 239–247.
35. I. Filella and J. Peñuelas, Daily, weekly and seasonal relationships among VOCs, NO_x and O₃ in a semi-urban area near Barcelona, J. Atmos. Chem., 2006, 54(2), 189–201.
36. C. Abadejo and L. Mariné, SOS: El Montseny s'ofega, Diari Ara, 2025.
37. J. Terradas, Holm Oak and Holm Oak Forests: An Introduction, 1999, pp. 3–14.
38. A. M. Yáñez-Serrano, A. Bach, D. Bartolomé-Català, V. Matthaios, R. Seco and J. Llusia, et al., Dynamics of volatile organic compounds in a western Mediterranean oak forest, Atmos. Environ., 2021, 257, 118447.
39. R. Seco, J. Peñuelas, I. Filella, J. Llusià, R. Molowny-Horas and S. Schallhart, et al., Contrasting winter and summer VOC mixing ratios at a forest site in the Western Mediterranean Basin: the effect of local biogenic emissions, Atmos. Chem. Phys., 2011, 11(24), 13161–13179.
40. R. Seco, J. Peñuelas, I. Filella, J. Llusia, S. Schallhart and A. Metzger, et al., Volatile organic compounds in the western Mediterranean basin: urban and rural winter measurements during the DAURE campaign, Atmos. Chem. Phys., 2013, 13(8), 4291–4306.
41. M. In 't Veld, A. Alastuey, M. Pandolfi, F. Amato, N. Pérez and C. Reche, et al., Compositional changes of PM2.5 in NE Spain during 2009–2018: a trend analysis of the chemical composition and source apportionment, Sci. Total Environ., 2021, 795, 148728.
42. A. Ripoll, M. C. Minguillón, J. Pey, N. Pérez, X. Querol and A. Alastuey, Joint analysis of continental and regional background environments in the western Mediterranean: PM1 and PM10 concentrations and composition, Atmos. Chem. Phys., 2015, 15(2), 1129–1145.
43. M. Cusack, A. Alastuey, N. Pérez, J. Pey and X. Querol, Trends of particulate matter (PM 2.5) and chemical composition at a regional background site in the Western Mediterranean over the last nine years (2002-2010), Atmos. Chem. Phys., 2012, 12(18), 8341–8357.
44. C. Selman and V. Misra, J. Geophys. Res., 2014,(3), 180–198.
45. N. Pérez, J. Pey, S. Castillo, M. Viana, A. Alastuey and X. Querol, Interpretation of the variability of levels of regional background aerosols in the Western Mediterranean, Sci. Total Environ., 2008, 407(1), 527–540.
46. M. Pandolfi, A. Alastuey, N. Pérez, C. Reche, I. Castro and V. Shatalov, et al., Trends analysis of PM source contributions and chemical tracers in NE Spain during 2004-2014: a multi-exponential approach, Atmos. Chem. Phys., 2016, 16(18), 11787–11805.
47. J. Pey, N. Perez, S. Castillo, M. Viana, T. Moreno, M. Pandolfi, J. M. López-Sebastián, A. Alastuey and X. Querol, Geochemistry of regional background aerosols in the Western Mediterranean, Atmos. Res., 2009, 94, 422–435.
48. M. In 't Veld, R. Seco, C. Reche, N. Pérez, A. Alastuey and M. Portillo-Estrada, et al., Identification of volatile organic compounds and their sources driving ozone and secondary organic aerosol formation in NE Spain, Sci. Total Environ., 2024, 906, 167159.
49. M. Ealo, A. Alastuey, N. Pérez, A. Ripoll, X. Querol and M. Pandolfi, Impact of aerosol particle sources on optical properties in urban, regional and remote areas in the north-western Mediterranean, Atmos. Chem. Phys., 2018, 18(2), 1149–1169.
50. M. Pandolfi, G. Martucci, X. Querol, A. Alastuey, F. Wilsenack and S. Frey, et al., Continuous atmospheric boundary layer observations in the coastal urban area of Barcelona during SAPUSS, Atmos. Chem. Phys., 2013, 13(9), 4983–4996.
51. A. Hansel, A. Jordan, C. Warneke, R. Holzinger, A. Wisthaler and W. Lindinger, Proton-transfer-reaction mass spectrometry (PTR-MS): on-line monitoring of volatile organic compounds at volume mixing ratios of a few pptv, Plasma Sources Sci. Technol., 1999, 8(2), 332–336.
52. L. Drinovec, G. Močnik, P. Zotter, A. S. H. Prévôt, C. Ruckstuhl and E. Coz, et al., The “dual-spot” aethalometer: an improved measurement of aerosol black carbon with real-time loading compensation, Atmos. Meas. Tech., 2015, 8(5), 1965–1979.
53. H. Hersbach, B. Bell, P. Berrisford, G. Biavati, A. Horányi, J. Muñoz Sabater, et al., ERA5 hourly data on single levels from 1940 to present, Copernicus Climate Change Service (C3S) Climate Data Store (CDS), 2023, https://cds.climate.copernicus.eu/datasets/reanalysis-era5-single-levels?tab=download.
54. P. Georgiades, T. Economou, Y. Proestos, J. Araya and J. Lelieveld, The Application of Machine Learning Algorithms to the Global Forecast of Temperature-Humidity Index with High Temporal Resolution, 2024.
55. A. Liaw and M. Wiener, Classification and Regression by randomForest, 2014.
56. M. Kuhn, Building predictive models in R using the caret package, J. Stat. Software, 2008, 28(5), 1–26.
57. A. Tani, Fragmentation and Reaction Rate Constants of Terpenoids Determined by Proton Transfer Reaction-Mass Spectrometry, Environ. Control Biol., 2013, 51(1), 23–29.
58. A. M. Yáñez-Serrano, I. Filella, J. LLusià, A. Gargallo-Garriga, V. Granda and E. Bourtsoukidis, et al., GLOVOCS - master compound assignment guide for proton transfer reaction mass spectrometry users, Atmos. Environ., 2021, 244, 117929.
59. J. L. Ambrose, K. Haase, R. S. Russo, Y. Zhou, M. L. White and E. K. Frinak, et al., A comparison of GC-FID and PTR-MS toluene measurements in ambient air under conditions of enhanced monoterpene loading, Atmos. Meas. Tech., 2010, 3(4), 959–980.
60. P. K. Misztal, M. R. Heal, E. Nemitz and J. N. Cape, Development of PTR-MS selectivity for structural isomers: monoterpenes as a case study, Int. J. Mass Spectrom., 2012, 310, 10–19.
61. R. Holzinger, L. Sandoval-Soto, S. Rottenberger, P. J. Crutzen and J. Kesselmeier, Emissions of volatile organic compounds from Quercus ilex L. measured by Proton Transfer Reaction Mass Spectrometry under different environmental conditions, J. Geophys. Res. Atmos., 2000, 105(D16), 20573–20579.
62. L. Fasbender, A. M. Yáñez-Serrano, J. Kreuzwieser, D. Dubbert and C. Werner, Real-time carbon allocation into biogenic volatile organic compounds (BVOCs) and respiratory carbon dioxide (CO₂) traced by PTR-TOF-MS, ¹³CO₂ laser spectroscopy and ¹³C-pyruvate labelling, PLoS One, 2018, 13(9), e0204398.
63. S. Owen, P. Harley, A. Guenther and C. Hewitt, Light dependency of VOC emissions from selected Mediterranean plant species, Atmos. Environ., 2002, 36(19), 3147–3159.
64. H. D. Castelyn, J. J. Appelgryn, M. S. Mafa, Z. A. Pretorius and B. Visser, Volatiles emitted by leaf rust infected wheat induce a defence response in exposed uninfected wheat seedlings, Australas. Plant Pathol., 2015, 44(2), 245–254.
65. M. L. White, R. S. Russo, Y. Zhou, J. L. Ambrose, K. Haase and E. K. Frinak, et al., Are biogenic emissions a significant source of summertime atmospheric toluene in the rural Northeastern United States?, Atmos. Chem. Phys., 2009, 9(1), 81–92.
66. A. Mellouki, T. J. Wallington and J. Chen, Atmospheric Chemistry of Oxygenated Volatile Organic Compounds: Impacts on Air Quality and Climate, Chem. Rev., 2015, 115(10), 3984–4014.
67. R. B. Stull, Introduction Boundary Layer-Meteorology, 1988, p. 671.
68. S. Kumar Mehta, M. Venkat Ratnam, S. V. Sunilkumar, D. Narayana Rao and B. V. Krishna Murthy, Diurnal variability of the atmospheric boundary layer height over a tropical station in the Indian monsoon region, Atmos. Chem. Phys., 2017, 17(1), 531–549.
69. S. Han, H. Bian, X. Tie, Y. Xie, M. Sun and A. Liu, Impact of nocturnal planetary boundary layer on urban air pollutants: measurements from a 250-m tower over Tianjin, China, J. Hazard. Mater., 2009, 162(1), 264–269.
70. C. Warneke, S. A. Mckeen, J. A. D. Gouw, P. D. Goldan, W. C. Kuster and J. S. Holloway, et al., Determination of urban volatile organic compound emission ratios and comparison with an emissions database, J. Geophys. Res. Atmos., 2007, 112, D10S47.
71. A. Gelencsér, K. Siszler and J. Hlavay, Toluene-benzene concentration ratio as a tool for characterizing the distance from vehicular emission sources, Environ. Sci. Technol., 1997, 31(10), 2869–2872.
72. S. Han, Z. Li, Y. Xiao, W. Wang, J. Li and X. Guo, et al., NO_x regime-dependent effects of biogenic emissions on toluene degradation and product formation, J. Environ. Sci., 2025 DOI:10.1016/j.jes.2025.10.008.
73. R. Atkinson and J. Arey, Atmospheric Degradation of Volatile Organic Compounds Atmospheric Degradation of Volatile Organic Compounds, Chem. Rev., 2003, 103, 4605–4638.
74. J. Llusia, J. Peñuelas, R. Seco and I. Filella, Seasonal changes in the daily emission rates of terpenes by Quercus ilex and the atmospheric concentrations of terpenes in the natural park of Montseny, NE Spain, J. Atmos. Chem., 2012, 69(3), 215–230.
75. R. Seco, J. Peñuelas and I. Filella, Short-chain oxygenated VOCs: emission and uptake by plants and atmospheric sources, sinks, and concentrations, Atmos. Environ., 2007, 41(12), 2477–2499.
76. Y. L. Dorokhov, E. V. Sheshukova and T. V. Komarova, Methanol in plant life, Front. Plant Sci., 2018, 871, 1–6.
77. K. Hüve, M. M. Christ, E. Kleist, R. Uerlings, Ü. Niinemets and A. Walter, et al., Simultaneous growth and emission measurements demonstrate an interactive control of methanol release by leaf expansion and stomata, J. Exp. Bot., 2007, 58(7), 1783–1793.
78. N. Dudareva, A. Klempien, J. K. Muhlemann and I. Kaplan, Biosynthesis, function and metabolic engineering of plant volatile organic compounds, New Phytol., 2013, 198(1), 16–32.
79. R. Atkinson and J. Arey, Gas-phase tropospheric chemistry of biogenic volatile organic compounds: a review, Atmos. Environ., 2003, 37(suppl. 2), S197–S219.
80. H. Insam and M. S. A. Seewald, Volatile organic compounds (VOCs) in soils, Biol. Fertil. Soils, 2010, 46(3), 199–213.
81. J. Peñuelas, D. Asensio, D. Tholl, K. Wenke, M. Rosenkranz and B. Piechulla, et al., Biogenic volatile emissions from the soil, Plant Cell Environ., 2014, 37(8), 1866–1891.
82. J. Peñuelas and J. Llusià, BVOCs: plant defense against climate warming?, Trends Plant Sci., 2003, 8(3), 105–109.

---