Contrasting organic aerosol molecular composition between the urban and agricultural environment of the Po Valley

The molecular composition of organic aerosol in the Po Valley remains largely unexplored, despite contributing approximately half of the fine aerosol mass. Molecularly-resolved analysis of the organic fraction is essential for understanding the sources and formation processes of organic aerosol in detail. Here, we investigated one year of PM_2.5 filter samples from a rural agricultural and an urban environment. We used liquid chromatography and high-resolution mass spectrometry with subsequent non-target analysis of 250 samples. Time-series analysis allowed for the grouping of detected organic compounds using a clustering algorithm, enabling a tentative source attribution. The most populated clusters consist mainly of CHO- and CHOS-containing compounds, attributed to oxidation products of biogenic emissions. They account 31 % and 26 % of the total intensity in the urban and agricultural sites, respectively, and peak during summer. Combustion-related clusters, enriched in nitrogen-containing compounds, contribute approximately 35 % of the total intensity at both sites. A fraction of these molecules are highly conjugated compounds that appear during winter as primary contributors to light-absorbing organic aerosol. Furthermore, we identified site-specific compound clusters, either at the urban or at the agricultural site. At the latter, we found pesticides strongly influence the overall molecular composition, peaking in May when PM_2.5 concentration is at its lowest level. This suggests potential toxicological effects despite apparent good air quality conditions. Our results represent the first molecular-level characterization of organic aerosol in the Po Valley, emphasizing the need to assess its composition for a better understanding of its environmental and health impacts.

1 Introduction

The chemical characterization of atmospheric aerosols is essential for developing effective mitigation strategies targeting the most critical sources of air pollution. It is well-established that aerosols affect human health, particularly among the most vulnerable groups of the population (Fann et al., 2018; Kioumourtzoglou et al., 2015; Mikati et al., 2018; Bell and Ebisu, 2012). High exposure of aerosols causes respiratory apparatus cancer (Chang et al., 2022), increases the risk for cardiovascular diseases (Wu et al., 2020; Liang et al., 2022; Chen et al., 2022; Peralta et al., 2022; Xie et al., 2018; He et al., 2022a), and decreases global cognitive performance (Weuve, 2012; Saenz et al., 2018; Ailshire et al., 2017; Harris et al., 2015) resulting in increased hospitalizations burdened on national health systems (Lee et al., 2019). Furthermore, aerosols interact with radiation, leading to a modification of Earth's energy balance. Ammonium nitrate and ammonium sulfate aerosols increase light scattering, which causes atmospheric cooling (Zhu et al., 2015; Adams et al., 2001). Drugé et al. (2019) estimated that ammonium nitrate aerosols cause a mean direct radiative effect at the top of the atmosphere of −1.4 W m⁻² over Europe, with a local maxima of −5 W m⁻² over the Po Valley, Italy, in all sky conditions. In contrast, black carbon generated by incomplete combustion is a strong light absorber and has a positive radiative effects (+0.0–4.0 W m⁻² over Europe, Nordmann et al., 2014), locally increasing air temperature (Ramanathan and Carmichael, 2008). In addition, other light-absorbing aerosols (LAAs) are organic aerosols (OA) formed by compounds rich in chromophores (Zeng et al., 2020), which strongly absorb radiation in the ultraviolet and visible region. These compounds are collectively referred to as “brown carbon” (BrC) due to the color that aerosol samples rich of this fraction assume. Ferrero et al. (2018) estimated BrC heating rate in the Po Valley accounts for 12.50 ± 0.06 % of the total, leading to a heating of 0.20 ± 0.01 K d⁻¹ in cirrostratus cloud conditions and 0.020 ±  10⁻³ K d⁻¹ under conditions of complete cloud cover (8 oktas) (Ferrero et al., 2021b). Molecular-level characterization of OA is becoming increasingly important to infer aerosol physicochemical properties, toxicological effects, and their sources (e.g. biogenic vs. anthropogenic, primary vs. secondary OA). This is undoubtedly challenging due to the number of compounds involved. Recently, Thoma et al. (2025) have investigated the OA fraction at the Taunus Observatory, a rural background station located at about 800 m a.s.l. in central Germany. By means of a multivariate analysis on an extensive dataset generated by an ultra-high-resolution mass spectrometer, the authors observed that even at a rural site, OA system is composed of several thousands of compounds affected by seasonality and short- and long-range transport, with biogenic secondary OA (BSOA) peaking during summertime and representing about 70 % of the total number of compounds, and 30 % attributed to anthropogenic secondary OA (ASOA), discerned because of the presence of combustion tracers, such as nitrophenols and aromatic species. The source attribution is even more important for air pollution hot-spot regions. Because of this, regions like the North China Plain (e.g., Zhang et al., 2023b; Zhu et al., 2021; Li et al., 2019) and the Indo Gangetic Plain (e.g., Panda et al., 2025; Gupta et al., 2022) have been extensively studied. For example, Kumar et al. (2022) reported that about 60 % of OA mass has a primary origin in New Delhi, due to hydrocarbon-like OA (HOA) and primary biomass burning OA (BBOA). Secondary organic aerosols (SOAs) were instead attributed to aromatics, BSOA, aged BBOA and mixed-urban SOA. Kapoor et al. (2024) pointed out the importance of LAAs, estimating that daily average BrC absorption contributes between 18 % and 42 % to near-ultraviolet absorption over India. Another well-studied pollution hot-spot is the Po Valley, in the Northern part of Italy. A number of papers (e.g., Neuberger et al., 2025; Bigi and Ghermandi, 2016; Ferrero et al., 2014, 2012) reported that this area is dramatically affected by its orography, which often prevents incoming air masses coming from North, West, and East able to disperse and abate the anthropogenic emission of the densely-populated valley in the lower layers of the troposphere. This leads to the accumulation of primary emitted compounds, especially during the cold season, when vertical mixing is limited to a shallow layer (Ferrero et al., 2011). In addition, Colombi et al. (2024) have shown that aerosol homogeneously affects the basin with very similar values in both concentration and chemical composition. The authors analyzed eleven years of particulate matter (PM) samples collected in an urban (Milan, the most populated city in the Po Valley) and rural (Schivenoglia, in a low-intensity agricultural region, about 150 km South-East of Milan) site, finding that both PM and secondary inorganic aerosols constituents show very similar concentrations along the year. Many studies have been published on source apportionment in the Po Valley, consistently highlighting the dominant role of anthropogenic emissions. For example, Perrone et al. (2012) performed a source apportionment on PM_2.5 samples collected at urban and rural sites over a three-year period (2006–2009), identifying the local contribution of traffic (Tr), the widespread influence of biomass burning (BB), and the importance of secondary aerosols (both organic and inorganic) particularly during spring and summer. Similarly, Bernardoni et al. (2013) used ¹⁴C analysis on PM₁₀ samples collected in Milan during wintertime to distinguish organic carbon (OC) and elemental carbon (EC) from fossil and modern carbon sources, estimating that OC_fossil accounts for about 30.8 % of total carbon, while primary wood burning contributes 14.8 %, with an additional 37.5 % attributed to other modern sources. Source apportionment of the OA fraction in this region has been extensively performed over the past decade using aerosol mass spectrometry (AMS). Several authors (Paglione et al., 2020; Saarikoski et al., 2012; Dall’Osto et al., 2015), conducted field campaigns at both urban and rural sites. Their results consistently showed that biomass burning organic aerosol (BBOA) significantly contributes to OA mass, especially in winter, through both primary and secondary components. Regarding the primary fraction, Decesari et al. (2014) highlighted the importance of hydrocarbon-like OA (HOA) and a cooking-related component, whereas Daellenbach et al. (2023) estimated that primary BBOA accounts for 37 % of the OA mass. Concerning the secondary fraction, SOA was found to dominate summertime, especially due to low-oxygenated compounds (Daellenbach et al., 2023). Sandrini et al. (2016) demonstrated that photochemical oxidation of gaseous precursors is a key mechanism for SOA formation, particularly under stagnant conditions. In addition, nighttime increases in relative humidity promote the formation of secondary water-soluble organic carbon, further enhancing the complexity of OA composition (Saarikoski et al., 2012; Decesari et al., 2014). These processes lead to the formation of light-absorbing “brown carbon” (BrC), often associated with the droplet mode (400–700 nm), enriched in aged OA and nitrate (Costabile et al., 2017). In this regards, Kiendler‐Scharr et al. (2016) showed that, in the southern Po Valley, the fraction of organic nitrates is comparable to that of inorganic nitrate and among the highest observed in Europe. Aqueous-phase reactions have been identified as a key pathway for BrC formation (Gilardoni et al., 2016), while a significant fraction of BrC has been found to be water-insoluble and associated with primary emissions from wood combustion (Gilardoni et al., 2020). Overall, the Po Valley atmosphere is strongly impacted by anthropogenic emissions, with more pronounced effects during winter, when meteorological conditions and orography favor air stagnation. Increasing dry periods throughout the year (Baronetti et al., 2020; Monteleone and Borzí, 2024) further enhance the accumulation and spatial homogenization of secondary species (Colombi et al., 2024). Consequently, source apportionment studies have reported relatively low spatial variability across the region, identifying traffic-related fossil emissions as a persistent source throughout the year, secondary aerosols as dominant in summer, and biomass burning as prevalent in winter. Despite these extensive efforts, a detailed molecular-level characterization of OA remains limited. Although AMS-based techniques are powerful tools for apportioning non-refractory OA sources via positive matrix factorization, their hard ionization approach limits the identification of specific molecular species. As a result, a comprehensive understanding of OA composition at the molecular scale is still lacking in this major air pollution hot-spot. Such knowledge is essential to improve the representation of OA in atmospheric models, particularly with respect to photochemical processes and parameterizations (Basla et al., 2022). Promising information to infer source attribution is provided by OA molecular characterization by ultra-high-resolution mass spectrometry (UHRMS) after soft ionization methods. Due to the high mass resolution and accuracy it is possible to derive molecular formulae from the exact mass. In addition, the chromatographic separation at the front of the mass spectrometer helps to resolve isomers. In this regard, many authors have already attempted to describe an aerosol system coupling ultra-high performance liquid chromatography (UHPLC) to soft-ionization UHRMS. For example, Mazzi et al. (2024) recently used a high-performance anion-exchange chromatograph coupled with a triple quadrupole mass spectrometer with an electrospray ion (ESI) source to detect and quantify plant protection agents in Mestre, near Venice, on the eastern side of the Po basin. Among the trace compounds found, the authors estimated glyphosate concentration of 0.05 ng cm⁻³ in PM₁₀. Many authors have begun to take advantage of the high resolution and accuracy of mass spectrometers, such as Orbitrap, to explore the complexity of OA without targeting specific molecules (e.g., Hildmann and Hoffmann, 2024; Song et al., 2023; Divisekara et al., 2023). Non-targeted analysis enables the characterization of the overall properties of OA and, when applied to time-series of samples and supported by multivariate analysis, has the potential to identify relationships between groups of compounds, molecular markers, and their sources (e.g., Wang et al., 2024). These relationships can be further constrained using source-specific analytical standards, external datasets, or other endpoints, as recently demonstrated by Breuninger et al. (2026) for the oxidative potential of PM. As noted by Tang et al. (2026), and originally proposed by Thoma et al. (2022), the concept of “Aerosolomics” advocates the integration of controlled experiments to elucidate chemical reaction pathways with atmospheric observations, time-series analysis, and correlations with known markers, supported by multivariate analysis, to identify the sources of atmospheric molecular species. In this context, the present work provides a substantial dataset of compounds that contributes to the Aerosolomics database, focusing on the Po Valley, a region for which, to our knowledge, these mentioned analytical techniques and data analysis approaches have not previously been applied. In this respect, Duarte et al. (2023) pointed out that, being the OA one of the most complex environmental matrices known, targeted and untargeted profiling of OA are both required to extract chemical information from the voluminous amount of data generated using newly available analytical approaches such as the one presented in this study. In this work, we applied the aforementioned techniques to two sets of fine atmospheric particulate matter samples (PM_2.5) collected in Milan (MI) and Schivenoglia (SKI), both within the Po Valley, during 2021. Although OC concentrations are frequently comparable between the agricultural (SKI) and urban (MI) monitoring sites (90° percentile of absolute difference: 2.0 µg m⁻³), the molecular characterization of OA reveals marked qualitative differences in molecular aerosol composition, resulting in different light-absorbance spectral behavior too. Furthermore, the implementation of a compound-level clustering analysis allowed for the source attribution of potentially light-absorbing BrC constituents, elucidating distinct source profiles between the two environments.

2 Method

2.1 Sampling

PM_2.5 was collected at two different locations, Milan-Pascal (45°28^′42.0348° N, 9°13^′54.03° E) and Schivenoglia (45°1^′0.336° N, 11°4^′34.0032° E), both belonging to the ARPA Lombardia (the regional environmental protection authority) air quality network. The first one (hereinafter “MI”) is an urban background site, located in a small park next to a low-traffic street of the University of Milan campus. Milan is the most populated city in the Northern part of Italy and is located in the center of the Po Valley. The second site (hereinafter “SKI”) is a rural background site about 150 km far from MI, surrounded by agricultural fields and few animal husbandries. A map of Northern Italy showing the sampling sites is provided in Fig. S1 in the Supplement. At both sites, low volume samplers (SKYpost PM, TCR-Tecora, or Lifetek PMS, Megasystem, 16.67 L min⁻¹) were used to collect the PM_2.5 samples on not pre-treated quartz fiber filters (Pall TissueQuartz, Ø = 47 mm) for 24 h since 2014. On a weekly basis, the filters were collected and stored in the darkness at low temperature (about 4 °C). For this work, we chose to investigate samples collected in 2021. Therefore, a total of 132 samples for MI site and 117 samples for SKI site were used for the chemical and molecular characterization. The covered period ranges from the end of January 2021 to the end of the same year. Moreover, additional data, such as meteorological parameters (air temperature, global radiation, relative humidity) and gaseous concentration of NO_x (NO + NO₂), benzene and toluene, were provided by ARPA Lombardia.

2.2 Chemical analysis

Water-soluble inorganic compounds and levoglucosan were quantified by high-performance anion-exchange chromatography with pulsed amperometric detection by extracting 1.5 cm² quartz filter punches in ultrapure water (10 mL volume) for 20 min in an ultrasonic bath. After filtration (Nylon or PTFE Syringe Filter, pore size 0.45 µm), the resulting solution was injected in the analytical system for the determination of anions (Metrohm 930), cations, and levoglucosan (Metrohm 881). A second punch was used to quantify EC and OC through thermo-optical analysis (Sunset Laboratory Inc., Tigard, OR, USA) based on the EUSAAR2 protocol. Both analytical procedures are also described in Colombi et al. (2024). For the organic phase investigation, a 12 mm diameter round punch was used for OA extraction in 180 µL of a 10 % acetonitrile (ACN) solution (Merck KGaA and Optima LC/MS Grade, Thermo Fisher Scientific Inc. and Milli-Q water Reference A+). The extraction was carried out using an orbital shaker (300 rpm) for 20 min. To the extracted and filtered (non-sterile PTFE Syringe Filter, Thermo Fisher Scientific Inc.) solution of 80 µL, an internal standard (10 % of the final extracted volume) consisting of isotopically labeled benzoic acid (${\mathrm{C}}_{\mathrm{6}}{\mathrm{H}}_{\mathrm{5}}^{\mathrm{13}}{\mathrm{CO}}_{\mathrm{2}}\mathrm{H}$, 99 at. % ¹³C, SigmaAldrich) and caffeine (¹³C₃C₅H₁₀N₄O₂, 99 at. % ¹³C SigmaAldrich) was added. The analytes were separated by injecting 10 µL of each sample into a UHPLC with a reverse phase column (CORTECS T3, 120 Å, 2.7 µm, 3 mm × 150 mm, Waters, C-18) system. The mobile phase, ultrapure water and ACN, was maintained at a constant flow rate of 0.4 mL min⁻¹ through the C-18 column heated at a constant temperature (40 °C). The gradient program was set as follows: 0–1 min at 1 % ($\mathrm{v}/\mathrm{v}$) ACN, 1–16 min with a linear increase to 99 % ACN, 16–16.5 min at a constant 99 % ACN, 16.5–17.5 min with a linear decrease back to 1 % ACN, and 17.5–20 min at 1 % ACN to re-equilibrate the column. A heated electrospray ion source (HESI-II Probe, Thermo Fisher Scientific Inc., operating in both positive and negative modes) was used to ionize the analytes before entering the mass spectrometer (high-resolution hybrid quadrupole-Orbitrap mass spectrometer, Q-Exactive Focus, Thermo Fisher Scientific Inc.). Ions in 75–750 $m/z$ range were detected in both positive and negative ionization mode in full-scan MS with a resolution of 70 000 at $m/z$ 200. Data-dependent tandem mass spectrometry (ddMS²) with a resolution of 17 500 was used at higher-energy collision energies (15, $\mathrm{41}/\mathrm{30}$ (negative/positive mode), 50 eV). As quality assurance, the sample analysis was periodically supplemented with mixed standards to evaluate the performance of the analytical system. A previous evaluation was also performed in Thoma et al. (2022).

2.3 Light-absorbance measurements

A subset of sample extracts was used to assess the light-absorbance to compare samples, which were collected at the two sites on the same days and exhibit very similar OC concentrations and contribution to PM_2.5 (within a tolerance of 2 µg m⁻³ and 5 %, respectively). For this purpose, 10 µL of sample extracts were directly injected without chromatographic separation into a diode array detector (DAD, Thermo Scientific, Vanquish, VH-D10-A) with a liquid carrier (99 % ACN, 300 µL min⁻¹). The light absorbance was measured with 20 Hz collection rate over the wavelength range of 230–650 nm. Solvent blank with the same amount of the internal standard injected in the sample extracts was used to remove the absorbance background. Finally, to account for the absence of a quantification of the OC mass extracted, we calculated the Absorbance Angstrom Exponent (AAE) for three wavelength ranges, that are 230–250, 250–400, and 400–650 nm. To this end, a linear interpolation was performed between the logarithm of absorbance and that of wavelength. The absolute value of the resulting slope represents the light-spectral behavior of the extracts.

2.4 Data analysis

2.4.1 OA non-target analysis

The chromatograms and acquired spectra were processed using Compound Discoverer software (Thermo Fisher, version 3.3.2). We considered spectra with retention time in the range 0–16.5 min and with a total intensity threshold of 10${}^{\mathrm{4}}/{\mathrm{10}}^{\mathrm{5}}$ (negative/positive mode). Retention time alignment was performed using an adaptive curve with a maximum absolute shift of 0.2 min. Each peak signal was extracted with a mass tolerance of ±2 ppm only if it appeared at least for 5 scans. Base ions for compound detection was set to [M − H]⁻ and [M + H]⁺ for negative and positive mode, respectively. Molecular formulae were attributed constraining the number of elements (min: C#1, H#1; max: C#90 H#190 Br#3 Cl#4 N#4 O#20 P#1 S#3), the DBE range (0–40), and the hydrogen-to-carbon ratio range ($\mathrm{H}/\mathrm{C}$, 0.1–3.5). The acquired MS² spectra were compared with the Advanced Mass Spectral Database mzCloud (HighChem LLC, Slovakia) and the Aerosolomic database, described in Thoma et al. (2022). A match with these libraries allowed to increase the identification confidence from Level 4 (L4, unequivocal molecular formula) to L2 (probable structure), based on the classification described in Schymanski et al. (2014). Detailed workflows are given in Tables S1 and S2 for negative and positive mode, respectively.

2.4.2 OA molecular characterization and hierarchical cluster analysis

After blank correction, the dataset was reduced considering only the compounds with a signal intensity (SI) higher than 10⁴ and 10⁵ in negative and positive mode, respectively, and with a signal-to-blank and a peak rate higher than 5. Moreover, to analyze the chromatographically resolved features only, we considered the retention time range between 1.55 min (void time: 1.44 min) and 16.50 min. The resulting list of molecular formulae were compared with the PubChem database (Kim et al., 2024, https://pubchem.ncbi.nlm.nih.gov/, last access: June 2024) to remove only potentially non-existent or impossible formulae. Assigned molecular formulae were used to group the features into families based on the elemental composition: CHO, CHNO, CHOS, CHNOS, CHN, and CHOP. Features without an assigned formula or not falling into the aforementioned categories were grouped under the category “Others”. To compare SI among the samples, each individual signal was scaled by the air volume sampled through the filter. Finally, to account for the different instrumental response in terms of signal peak intensity due to the ionization modes, the SI were normalized (nSI) to the maximum value found for each polarity (see Eq. S1). The non-target screening approach requires the identification of descriptors to infer the properties or characteristics of detected molecules, which can aid in their characterization. For each formula, we calculated the average carbon oxidation state (OS_c) with the simplified method described in Kroll et al. (2011) (Eq.1), the double bond equivalent value (DBE, Eq. 2), that is the number of rings and double bonds in a molecule, the aromaticity equivalent (X_c), as proposed in Yassine et al. (2014) (Eq. 3). This last descriptor was validated to provide information concerning the presence of aromatic structure in a molecule (X_c ≥ 2.5000), while aliphatic structures show lower values of the mentioned threshold. In this work, to infer temporal averages, we used the peak abundance-weighted average parameter approach, as already shown in Wang et al. (2021).

$$\begin{array}{}\text{(1)}& {\displaystyle }{\mathrm{OS}}_{\mathrm{c}}=\mathrm{2}\times \mathrm{\#}\mathrm{O}/\mathrm{\#}\mathrm{C}-\mathrm{\#}\mathrm{H}/\mathrm{\#}\mathrm{C}\text{(2)}& {\displaystyle }\mathrm{DBE}=\mathit{\#}\mathrm{C}-\mathit{\#}\mathrm{H}/\mathrm{2}+\mathit{\#}\mathrm{N}/\mathrm{2}+\mathrm{1}\text{(3)}& {\displaystyle }{X}_{\mathrm{c}}={\displaystyle \frac{\mathrm{3}\times [\mathrm{DBE}-(p\times \mathit{\#}\mathrm{O}+q\times \mathit{\#}\mathrm{S}\left)\right]-\mathrm{2}}{\mathrm{DBE}-(p\times \mathit{\#}\mathrm{O}+q\times \mathit{\#}\mathrm{S})}}\end{array}$$

In Eq. (3), p and q represents the fraction of oxygen and sulfur atoms involved in the π-bond structures. Since their values depend on the category of compounds, we used the same approach proposed in Kourtchev et al. (2016),Wang et al. (2017), and Tong et al. (2021), where $p=q=\mathrm{0.5}$ for compounds detected in HESI(−), while $p=q=\mathrm{1}$ for formulae detected in HESI(+). In addition, the value of (p× #O + q × #S) is rounded down to the lower integer whether #O or #S are odd for compounds detected in negative mode according to Yassine et al. (2014). It is worth stating here as well that the analytical approach used in this work is not suitable for the quantification of each compound. Although the SI is proportional to the concentration in the analyzed sample, soft ionization techniques have different ionization efficiencies for each analyte. Theoretically, it would be possible to measure the ionization efficiency of each compound; however, the enormous number of molecules contained in an atmospheric aerosol sample makes this process practically impossible. Moreover, many compounds present in the atmosphere are completely unknown and, therefore, no analytical standard is available. This issue led some authors (e.g., Ma et al., 2022; Divisekara et al., 2023) to use surrogate compounds to represent a broader group of molecules and to assume the same concentration-peak area intensity response for these as for the measured compound. Although, to our knowledge, this is the best approach currently available, we believe that the number of surrogates required would be too large for the purposes of this work. For this reason, we have limited our discussion to considering peak area intensities as a proxy for atmospheric concentration, i.e. to compare their relative importance in the samples. A hierarchical clustering algorithm was chosen to reduce the complexity of the system (Qi et al., 2019; Priestley et al., 2021; Thoma et al., 2022). We used the software MATLAB (The MathWorks, vers. 2024a) independently on the two datasets after removing all the features with an occurrence lower than 10 %. We standardized (through centering) each feature along its time-series and used the Euclidean distance as a metric to calculate the similarity among the compounds and the Ward linkage method to create the tree. Finally, the number of clusters for each site was chosen based on the similarity metric, the resulting fingerprint of each cluster, its interpretability using available markers in the literature, and its correlation with external data, such as EC, soluble potassium, levoglucosan, and NO₂.

3 Results

3.1 Seasonality of organic aerosol molecular composition

The calculation of the normalized total signal intensity (nTSI, Eq. S2), i.e. the sum of the nSI of individual compounds in each sample, allowed us to compare the seasonal variability at each site. Both time series sites exhibit peak values during winter and summer (Fig. 1). A similar trend is also visible in the OC concentration quantified in the PM_2.5 samples from both sites, although the increase during the summer period, mainly due to aliphatic CHO and CHOS compounds, is not fully reflected in the OC values. This leads to a poor correlation (R² = 0.280 and 0.466 for MI and SKI, respectively) between nTSI and OC on an individual samples basis. Highly oxygenated compounds dominate during summer, with an average OS_c of 0.91 ± 0.06 for MI and 0.89 ± 0.04 for SKI in July. Figure 1 indicates that nTSI has a stronger contribution in aromatic compounds, especially belonging to CHNO and CHO families, during winter. In this regards, both DBE and X_c index show minimum values during summer (2.54 ± 0.14 and 0.37 ± 0.09, respectively in July) and maxima in winter (3.40 ± 0.34 and 0.87 ± 0.15, in February). While no clear temporal trend is observed for P-containing compounds in Milan, two periods of notable contribution from this family to the nTSI are identified at the agricultural site. The annual pattern shows that this compounds contribution increases from March to early May and from late September to the beginning of December, peaking in samples collected in April and November.

Figure 1 Time series of the monthly average nTSI for the groups of compounds detected in the MI (top plot) and SKI (bottom plot) samples. Darker colors represent the nTSI of aromatic compounds (“−a”), distinguishing them from the lighter colors for aliphatic ones (“−n”) based on their X_c values.

Retrieving light-absorption information via a molecular-derived proxy, we use the $\mathrm{DBE}/\mathrm{\#}\mathrm{C}$ ratio to identify compounds of light-absorbing aerosols (LAAs), setting limits for this range between 0.5 (for polyenes) and 0.9 (for fullerene-like hydrocarbons) (e.g., Lin et al., 2018; Laskin et al., 2014). Using these criteria, we selected molecules from our dataset that likely contribute to BrC (Fig. 2). Organic LAAs show maxima nTSI values during winter (0.73 ± 0.36 at MI and 0.53 ± 0.24 at SKI). This is in agreement with Bluvshtein et al. (2017) and Budisulistiorini et al. (2017) describing that emission from wildfires, peat fires, agricultural residue burning, and residential combustion substantially contribute to BrC. While some of these sources may remain active during summertime, photobleaching and ventilation of the Po Valley are more efficient during summertime (Chen et al., 2021). Consistently, we observe lower values of the BrC-related molecules in summer (0.27 ± 0.11 at MI and 0.16 ± 0.04 at SKI). The higher values in nTSI observed at the urban site in the warm season indicate either more active emission sources or air-quality conditions that more efficiently promote chromophores formation. Indeed, the plots in Fig. 2 show that the urban site (left) exhibits a larger number of BrC-related compounds (+409) compared to the rural site (right). From a qualitative perspective, the seasonal BrC plots for the two sites (Fig. S2) confirm a different composition of LAAs. Seasonally averaged BrC compound intensities were aggregated by molecular-family and normalized to the total BrC intensity (Fig. S3), revealing that the relative contributions of these families differ between the two sites across seasons. During them (summer and winter) the urban site exhibits higher contributions due to the CHN family compared to the agricultural site (8.8 % vs. 3.7 % at MI and SKI, respectively, in summer; 17.1 % vs. 7.9 % in winter), whereas CHO family contributes to the higher nTSI values (Fig. 2) in respect to SKI during summer only. In contrast, CHNO species contribute more at the agricultural site in both seasons (23.9 % vs. 33.9 % in summer and 46.2 % vs. 51.5 % in winter, respectively). Contributions from the other families are marginal, accounting for less than 2 % at both sites.. Similarly to the temporal pattern shown in Fig. 1, at both sites the contribution of wintertime BrC is mainly due to N-containing compounds. At the same time, the nSI-weighted average shows that in the coldest months BrC has the lowest $\mathrm{O}/\mathrm{C}$ ratio, suggesting a low oxidation regime affecting BrC molecular composition (Fig. S4). On the contrary, the lower summertime peak in BrC nTSI is mainly ascribable to highly oxidized compounds ($\mathrm{O}/\mathrm{C}$ > 0.5) with a probable secondary origin (Gilardoni et al., 2020).

Figure 2 Fingerprints of potential BrC compounds at MI (left) and SKI (right) sites and their nTSI time series over the year 2021. The orange-shaded area delineates the region in which potential BrC compounds are located.

In addition to the seasonality of elemental-composition groups and the BrC-related molecules, we focus on few individual tracers to confirm our hypothesis concerning their seasonal patterns. For instance, the later-generation product of monoterpene oxidation, that is 3-methyl-1,2,3-butane-tricarboxylic acid (MBTCA), peaks in summer. This is consistent with the higher monoterpene emission and stronger oxidative conditions from May to September. On the opposite, the CHNO compound nitrocatechol, a well-known combustion product from BB (Fredrickson et al., 2022), shows the opposite trend. The time series of these compounds, whose identification was allowed by the comparison with analytical standards (identification level 1, L1, based on Schymanski et al., 2014), are shown in Fig. S5.

Based on the blank-corrected annual average SI, we ranked the compounds with the highest values for the most representative molecular families (that are CHO and CHNO), which alone represent 60 % of the nTSI. Within the CHO(−) family, the formulae C₄H₆O₅, C₈H₁₂O₆, and C₆H₈O₇ exhibit the highest average intensities. The first two compounds show a distinct seasonal trend, with intensities peaking in spring and reaching a maximum during the summer, followed by a decline during autumn. This pattern suggests that the conditions favorable for the emission or formation of these compounds are most prominent in summer. Together with the aforementioned MBTCA (L1, C₈H₁₂O₆), the compound with the formula C₄H₆O₅ and its pronounced increase during summer suggest a likely monoterpene-derived biogenic origin (L4, as suggested by Thoma et al., 2022, Chen et al., 2020, and the literature within). In contrast, the third-highest intensity compound, C₆H₈O₇, appears to be an anthropogenic SOA tracer (ASOA) derived from the oxidation of benzene (L4, Cheng et al., 2021). For the CHNO(−) molecular family the highest signals occur during wintertime. At both sites, nitrocatechol (L1, C₆H₅NO₄) ranks highest. Additionally, at the urban site, two isomers of C₇H₇NO₄ were detected, likely representing methyl-nitrocatechol and nitroguaiacol (L4, Al-Naiema and Stone, 2017). Finally, SKI samples reveal a significant signal for C₈H₆N₂O₂, which can be attributed to nitroindole (L4), known as a BrC-related compound from the oxidation of indole in presence of NO₂ (Jiang et al., 2024). The presence of this gaseous amine in the atmosphere has been attributed to many sources, such as vegetation, biomass burning, industrial activities and animal husbandry (literature within Xue et al., 2022). We observe a similar annual pattern for nitroindole at urban (MI) and rural-agricultural (SKI) sites. However, SKI exhibits higher signal intensities, indicating a closer proximity to the precursor sources, i.e. animal husbandry. Since the $\mathrm{DBE}/\mathrm{\#}\mathrm{C}$ of nitroindole is 0.875, it is also highlighted in Fig. 2 as the blue circle with coordinates (7; 8). In positive ionization mode, the highest intensities for CHO(+) compounds are measured in MI as C₁₁H₂₂O₅ and C₁₂H₂₆O₅, attributed to ASOA from alkane photooxidation (Zhang et al., 2014), and more specifically to n-dodecane SOA (L4, Li et al., 2021). The third highest signal at MI corresponds to the highest one at SKI, i.e. C₁₀H₂₂O₄, although no correspondence was found in the literature. At SKI site, C₈H₁₆O₂ and C₂₀H₂₆O₃, the second and third highest signals, are attributed to BBOA by Kong et al. (2021) and Smith et al. (2009). For the nitrogen-containing compounds (CHNO) detected in positive mode, the formulae corresponding to the $m/z$ values with the highest signals,are C₉H₁₁NO₂, C₁₂H₁₇NO and C₆H₁₅NO for MI site and C₉H₁₈N₂O, C₁₂H₂₃NO and C₁₂H₂₅NO for SKI, tentatively assigned to 1-[[(2R)-oxolan-2-yl]methyl]piperazine, a plant protection product (L4, ECHA CHEM database, 2024b), a compound found in agricultural residue burning (L4, Li et al., 2024) or to n-octyl pyrrolidone (L4, ECHA CHEM database, 2024a) also used as co-solvent in crop protection formulations, and N,N-dimethyl capramide, an additive in pesticide formulations (L4, PubChem, 2024). As these last mentioned compounds are anthropogenic chemicals, an unambiguous molecular identification can be achieved via target analysis.

3.2 Organic aerosol hierarchical clustering

In the following, we describe the hierarchical cluster analysis of the complete dataset. We conducted two independent cluster analyses for the two sites to ensure that the differences between the sites are maximized, enabling an unbiased assignment of a compound to a cluster regardless of its temporal behavior at the other site. It is important to clarify that the presence of a formula within a cluster does not necessarily indicate that the corresponding compound is emitted by the source after which the cluster is named. Rather, it may suggest a shared emission driver or reactive pathway that causes the compounds in the same cluster to exhibit a similar temporal pattern. However, many anthropogenic sources are known to be active year-round, with their influence often varying over time. In addition to the seasonality of emission strength, the meteorological variability can drive different reaction pathways. An increased oxidation of anthropogenic VOCs during the summer months might consequently lead to the assignment to a biogenic cluster. Even more difficult to assign to an anthropogenic cluster are SOA compounds that originate from the oxidation of volatile chemical products (VCP), such as limonene. On the other hand, such products may experience an increase in their concentration during wintertime due to a reduced mixing layer height. As such, the following discussion does not intend to definitively attribute a source to each detected compound, but rather to provide insights into their temporal trends and which other compounds show similar behavior over the course of the year and assign them to a representative molecular or sector-specific cluster. Furthermore, the use of highly resolved chromatographic and spectrometric techniques yields a large number of isomers, which cannot be precisely identified in the absence of specific reference standards of characteristic fragmentation spectra. As a result, comparisons with data from literature may also be subject to this inherent uncertainty. At each site we examined the compounds that were detected in more than 10 % of the analyzed samples. This leads to investigate via hierarchical cluster analysis (HCA) 6487 and 5499 parameters for MI and SKI, respectively. We also used external data from meteorological factors (such as global radiation, daily temperature, and relative humidity), gaseous species concentrations (NO_x, O₃, benzene and toluene), ion chromatographic and thermo-optical transmittance data (i.e., EC and OC) to improve our interpretation of the results. The HCA revealed seven main clusters for the urban site (Fig. S6) and five for the rural site (Fig. S7). According to the temporal patterns, three clusters at both sites showed a clear increase in intensity during the colder season. The higher signals observed during winter may arise from multiple contributing factors, acting either independently or in combination. Seasonal variations in emissions may lead to the presence of sources that are active or more intense in winter than in other periods of the year. In addition, meteorological conditions typical of winter, such as reduced atmospheric mixing and dispersion, can enhance aerosol concentrations. Partitioning processes may also play a role, with gas–particle equilibrium shifting toward the particulate phase under colder conditions. Furthermore, the absence or low abundance of certain compounds in summer does not necessarily imply a lack of emissions; rather, enhanced photochemical degradation may lead to their rapid transformation into other species. As a result, compounds that are prevalent in winter may be replaced in summer by secondary products that are only present at low levels during the colder season. In either case, this behavior suggests that the compounds in these clusters have a likely anthropogenic origin. Additionally, at MI, two clusters show an increase during the warmer season, whereas only one cluster exhibits this behavior at the rural site. The remaining clusters do not show a clear temporal pattern, which could be due to a sporadic rise or activation of the sources. The contribution of the clusters, as monthly average, is shown in Fig. 3. In the Supplement, we show the fingerprints of each cluster using the retention time versus the molecular mass of each compound, the Van Krevelen diagram, and the Kroll plot to infer the cluster's chemical signature. Finally, the molecular formulae with tentative identifications and cluster attribution are listed in Additional Data (Table A1).

Figure 3 Time series of the nTSI of the clusters identified at MI (upper plot) and SKI (lower plot) sites.

3.2.1 Wintertime clusters

The molecular fingerprints of the wintertime clusters (Figs. S8, S9, S10, S15, S16, and S17) are characterized by a significant contribution of N-containing compounds and low $\mathrm{O}/\mathrm{C}$ ratio values (median < 0.45). The region of the Van Krevelen diagram that is typically associated with aromatic compounds ($\mathrm{O}/\mathrm{C}$ between 0 and approximately 1 and $\mathrm{H}/\mathrm{C}$ between approximately 0.5 and 1) is densely populated. At the urban site, we identify a distinct cluster in which this region is predominantly occupied by CHO compounds (Fig. S8). Notably, the nTSI of this cluster strongly correlates with soluble potassium (R²= 0.705) and levoglucosan (R²= 0.734), both established tracers of biomass burning sources (e.g., Yu et al., 2018). This suggests that biomass burning is likely the dominant source for most compounds in this cluster. A few CHNO compounds, such as nitrophenol (C₆H₅NO₃, L1), isomers of methyl-nitrophenol (C₇H₇NO₃, L4), methyl-nitrocatechol or nitroguaiacol (C₇H₇NO₄, L4), vanillin (C₈H₈O₃, L4) and vanillic acid (C₈H₈O₄, L4), also fall within the aromatic region of the Van Krevelen diagram of this cluster. In the CHO class, we identified phthalic acid (L1) and its isomers. Phthalic acid is released into the atmosphere from several anthropogenic sources, such as combustion of plastics or oxidation of traffic exhaust PAHs, such as naphtalene (e.g., Lui et al., 2023; Al-Naiema et al., 2020; Al-Naiema and Stone, 2017; Deshmukh et al., 2016; Kawamura and Kaplan, 1987). The presence of biomass burning markers and the correlation with the external tracers suggest that BB is the dominant source for the “MI-BBOA” cluster. Its nTSI has the highest coefficient of determination with OC concentrations (R²=0.813), suggesting also that BB is a relevant source for carbonaceous matter concentrations at MI site during wintertime. The nTSI pattern of the cluster called “MI-Other combustion” (Fig. S9) shows a good correlation with soluble potassium (R²= 0.718) and levoglucosan (R²= 0.794, populated by nitrocatechol (C₆H₅NO₄,L1), three isomers of methyl-nitrocatechol (C₇H₇NO₄, L4) and coniferyl alcohol (C₁₀H₁₂O₃, L1). The characteristic of this cluster is mainly the presence of a large number of CHN(+) compounds, such as C₇H₁₂N₂, C₈H₁₄N₂, C₉H₁₀N₂, C₁₂H₁₀N₂, which are formulae attributed to BBOA by Smith et al. (2009) and Fleming et al. (2018). Moreover, even with slightly weaker correlation, this cluster indicates an influence from traffic sources, as suggested by NO₂ (R²= 0.567) and benzene (R²= 0.570), leading to a cluster possibly affected by both sources. However, we also observe a weak correlation to EC_traffic (R²= 0.354), estimated by the difference between EC and EC_non-fossil (Daellenbach et al., 2023). This approach derives from the parametrization described in Zotter et al. (2014), which reported a EC_non-fossil–to–levoglucosan ratio of 0.87 for stations south of the Alps and assuming that EC is due to fossil (coal combustion and traffic exhaust) and non-fossil combustion (wood burning). On the other hand, these external tracers correlate better with a third wintertime cluster found at the urban site (R²= 0.648, 0.829 and 0.747 for EC_traffic, NO₂ and benzene, respectively), which fingerprint is shown in Fig. S10. The nTSI of this cluster shows increasing values during the weekdays and lower ones in the weekends. Many tracers attributed to traffic exhaust in the literature are found here: e.g. C₄H₆O₄ (succinic acid, L4, Lui et al., 2023), C₇H₁₂O₇ (L4, Thoma et al., 2022), and C₉H₈O₄ (methylphthalic acid, L4, Ikemori et al., 2021). In addition, N-containing features enrich the so-called “MI-TrOA” (Traffic OA) cluster, including aliphatic amines such as C₈H₁₈N, C₇H₁₇N, C₉H₂₁N (as reported in Cao et al., 2023). This cluster contains probable cyclic compounds, like C₁₀H₁₄N₂, which could be nicotine (L4) and nitro-aromatic compounds, C₆H₃N₃O₇ (picric acid, L4, Lin et al., 2015). More in general, it includes N-organic containing compounds (NOCs), some reported in the literature from traffic exhaust, like C₉H₁₉NO (N,N-dibutylformamide, Rogge et al., 1993, L4). These characteristics clearly distinguish this cluster from MI-BBOA and MI-Other Combustions clusters discussed earlier. In addition, we found some CHOP compounds, such as C₁₈H₃₉O₇P (tris(2-butoxyethyl)-phosphate, L4, Ungeheuer et al., 2021), C₁₈H₁₅OP (triphenylphosphine oxide, L4, Faiz et al., 2018), and C₈H₁₉O₅P (found in Lao et al., 2022), attributable to flame retardants. This class of compounds are often observed as non-combustion emission by vehicles too (Dong et al., 2025). Within the wintertime clusters, 96 out of the 239 compounds likely HOA are observed, based on their values of OS_c (between −2 and −1) and #C (between 18 and 30) (Kroll et al., 2011), distributed as 16 %, 30 % and 54 % among BBOA, Other Combustion and TrOA clusters, respectively. As an example, C₁₉H₂₆O₂ and C₂₀H₂₈O₃ are attributed to BBOA (Zhang et al., 2023a and Smith et al., 2009, L4), whereas C₁₈H₂₈O to vehicular emissions (Li et al., 2021, L4). At the rural site, a total of 82 HOA compounds are found, of which 46 appear in one of the three clusters with a winter temporal pattern. The cluster named “SKI-BBOA/TrOA” (Fig. S15) exhibits a strong correlation with both biomass burning (i.e. soluble potassium, R²= 0.674, levoglucosan, R²= 0.708) and traffic (i.e. benzene, R²= 0.769, EC_traffic, R²= 0.539, and NO₂, R²= 0.678) markers. Notable compounds include vanillic acid (C₈H₈O₈, L4, Simoneit, 2002), nitrophenol (L1), dodecanitrile (C₁₂H₂₃N, L4, Simoneit et al., 2003), nitrocatechol (L1), nitrosalicylic acid (C₇H₅NO₅, L4, He et al., 2022b), acetosyringone (C₁₀H₁₂O₄, L4, Li et al., 2020b and references therein), phthalic acid (L1), methylphthalic acid (L4), and 4-hydroxy-3-nitrobenzyl alcohol (C₇H₇NO₄, L4, Al-Naiema and Stone, 2017). Also, 33 HOA compounds fall in this cluster, such as C₂₀H₂₈O₃, and C₂₀H₂₆O₂, attributed by Smith et al. (2009) to BB sources. Similar to clusters identified in the MI dataset (Figs. S8 and S10), S-containing compounds were also detected in wintertime clusters. Among these, we found C₁₂H₂₀O₇S (L4, literature in Srivastava et al., 2022), C₈H₁₈O₄S (L4, Blair et al., 2017), both linked to SOA from traffic emissions or BSOA. These last include both polar compounds like C₅H₈O₇S (L4, SOA from isoprene and SO₂, Claeys and Maenhaut, 2021) or C₄H₆O₇S (L4, SOA from isoprene and SO₂, Lin et al., 2022), as well as larger compounds that were retained more effectively by the chromatographic column, such as C₁₀H₁₇NO₇S (L4, Surratt et al., 2008). The temporal pattern of these wintertime S-containing compounds suggests the importance of the co-presence of their biogenic precursors and high NOx concentrations. At the SKI site, wintertime S-containing compounds are equally distributed between SKI-BBOA/TrOA and the second cluster with a winter temporal pattern (Fig. S17). We observed a moderate correlation with inorganic nitrate (R²= 0.541) and ammonium (R²= 0.415) ion concentrations, known for their significant contribution to PM_2.5 during winter months in the Po Valley (Ferrero et al., 2019), EC (R²= 0.669, $R_{\text{EC-traffic}}^{\mathrm{2}}=$ 0.524) and benzene (R²= 0.515). On the other hand, the cluster fingerprint clearly shows a strong contribution of nitro-aromatics compounds, such as nitrophenol and nitrocatechol, indicating a probable combustion source (BB) contributing to OA with reaction products sensitive to the same air conditions for ammonium nitrate gas-to-particle partitioning, and OC concentrations (R²= 0.681). Because of the strong presence of nitro-aromatics, we labeled this cluster “SKI-Nitro Aromatics”. Within this cluster, 11 P-containing compounds significantly contribute to the CHOP family's nTSI, as shown in Fig. 1. Among these, C₆H₁₅O₄P, C₂H₇O₄P, and C₄H₁₁O₄P, show the highest nSI. The molecular formulae match with triethylphosphate, used as a flame retardant and as an additive in plant protection formulations, dimethylphosphate, and diethyl hydrogen phosphate, likely degradation products of other organophosphates (Liu et al., 2021). The third wintertime cluster detected at SKI (Fig. S16) contains many compounds assigned to “MI-Other combustion”. This SKI cluster primarily consist of NOCs (approximately 90 % of the features), mainly detected in positive ionization mode. This suggests a prevalence of reduced compounds, likely amine-related CHN(+). The X_c values below 2.5 for 65 % (median 2.3) of the features indicate that these compounds are mainly aliphatic. CHNO(+) compounds exhibit low $\mathrm{O}/\mathrm{C}$ (0.17 ± 0.04) and $\mathrm{O}/\mathrm{N}$ (1.2 ± 0.4) ratios, consistent with the presence of amino groups and O-containing substituents (average carbon number 8.11 ± 0.57). The peak abundance-weighted average of the nitrogen number in the CHN family is 1.98 ± 0.14. Based on these observations, we labeled this cluster as “SKI-Amines related”. The time series of this cluster reveals certain samples with high nTSI values for all the compounds it contains. This behavior is not observed for compounds in other clusters. The Kendrick mass defect analysis highlights homologous series such as C_nH_2nN₂ ($n=\mathrm{7},\mathrm{9},\mathrm{11}$), C_nH_2n−2N₂ (n= 5–13), C_nH_2nN₂ ($n=\mathrm{6},\mathrm{8},\mathrm{10},\mathrm{12}$). These series suggest probable methylation/demethylation processes in their genesis. However, previous studies have shown that the presence of reduced NOCs in the atmosphere can stem from a variety of sources. Murphy et al. (2007) reported that cooking, industrial sources, biomass combustion, and traffic are just a few of the possible sources, emphasizing that livestock activities are the predominant contributors of these compounds. In this context, it is worth to notice that other authors (such as Silva et al., 2008 and Malloy et al., 2009) have demonstrated that the reaction of aliphatic amines with nitrate radicals leads to the formation of imines in the aerosol phase, while Smith et al. (2021) have described the ability of imines to interact with light, suggesting their role as a component of secondary BrC.

3.2.2 Summertime clusters

At both sites we detect highly populated cluster peaking during summertime, with 1487 and 1623 compounds for MI and SKI, respectively, mostly detected in negative mode. The molecular fingerprints of these two clusters (Figs. S11 and S18) are very similar. Summertime clusters are mainly constituted by aliphatic CHO and CHOS compounds, with short carbon chains (peak-abundance weighted average, #C = 6.7 ± 0.7). Interestingly, the values of the averaged carbon number (weighted by the peak area) increase during the summertime more than the averaged oxygen number, leading to lower $\mathrm{O}/\mathrm{C}$ ratios. Because of this, the average OS_c decreases as well in the warmer season. Compared to the other summertime clusters, these clusters exhibit the highest $\mathrm{O}/\mathrm{C}$ ratio. CHO compounds have $\mathrm{H}/\mathrm{C}$ ratios ranging from 1 to 2 and $\mathrm{O}/\mathrm{C}$ ratios between 0.3 and 1.5. The high oxygen atom content (median = 6.0) observed in the S-containing compounds suggests that they are organosulfates. In this regard, many formulae match with ones we found in literature from isoprene oxidation products in presence of sulfuric acid, such as C₄H₈O₇S, C₄H₈O₆S (likely 2-methylglyceric acid organosulfate, L4, Claeys and Maenhaut, 2021), C₅H₁₂O₇S (SOA from 2-methyltetrol organosulfate, L4, Lin et al., 2022) and C₅H₁₀O₇S (L4, Nestorowicz et al., 2018), and from monoterpenes like C₈H₁₄O₇S, C₇H₁₂O₆S, C₁₀H₁₈O₇S (L4, literature within Chen et al., 2018; Jiang et al., 2024). The two clusters include MBTCA (L1), many isomers of C₈H₁₄O₄, which could be for example caric acid (L4, Gómez-González et al., 2012), SOA from limonene (C₇H₁₀O₅, L4), β-pinene (C₈H₁₄O₅) and Δ³-carene (C₆H₁₀O₂), also found by Thoma et al. (2022) using HESI(−). Based on these indicators, we conclude that the two clusters at MI and SKI represent summertime biogenic SOA (“BSOA”). Alongside the CHO and CHOS compounds, a considerable amount (that are 205 and 104 compounds, which represent 13.8 % and 6.5 % of the total compounds in BSOA cluster for MI and SKI, respectively), even with lower nTSI, are N-containing compounds detected in positive ion mode. Many of these were also listed in Zhang et al. (2023a), and attributed to wood burning (L4). Although the identification level is only 4 (unequivocal molecular formula), we hypothesize that these compounds are rising during summertime from oxidized BBOA. In this season, the oxidative capacity of the atmosphere is enhanced by higher solar radiation, which leads to higher concentrations of oxidants species both during daytime (mainly OH radicals and O₃) and nighttime (NO₃ radicals). In this regards, due to the 24 h time-resolution, a contribution of both daytime and nighttime oxidation is expected to alter the OA composition. The presence of BBOA during summer has been previously demonstrated (Ferrero et al., 2018, 2021a, b). Indeed, the moderate correlation between NOCs for MI-BSOA and daily average O₃ concentration (R²= 0.583) provides an insight of this. This hypothesis is further supported by the dominance of the CHNO(+) family contributing to the nTSI of the second summertime cluster identified at the urban site (Fig. S12). The whole group is mainly populated by compounds detected in positive ion mode and belonging to CHO and CHNO families (329 and 225 features, respectively). Both exhibit higher correlation with daily average O₃ (R²= 0.576 and 0.458, respectively) and global radiation (0.434 and 0.381) values in respect to the same families belonging to MI-BSOA (0.414 and 0.366 with O₃, 0.212 and 0.190 with the global radiation). While the cluster's $\mathrm{H}/\mathrm{C}$ ratio is close to the one of MI-BSOA (median 1.63 vs. 1.56, respectively), its $\mathrm{O}/\mathrm{C}$ ratio is significantly lower (0.37 vs. 0.60), resulting in a lower averaged OS_c (−0.9 vs. −0.3). In this cluster we also found many features attributed to biogenic sources, such as C₇H₁₀O₄, C₈H₁₁NO₇, C₁₀H₁₅NO₉ attributed to SOA from monoterpenes (L4, Finessi et al., 2014; Faxon et al., 2018; Pullinen et al., 2020), and falling in the SKI-BSOA cluster too. However, more than 40 % are features not detected at SKI, of which 70 % are detected in positive mode, like C₉H₂₀O₄ and C₁₃H₂₄O₄, which instead are found in Qi et al. (2019) in a urban site in Zurich and apportioned to cigarette smoke OA and cooking OA, respectively. Because of this strong mix of both biogenic and anthropogenic source assignments, we labeled this cluster as “MI-Summer B/ASOA”. Although terpenes are traditionally considered volatile organic compounds of predominantly biogenic origin, it is well established that fragranced consumer products, including cleaning agents, air fresheners, and personal care products, can contain and emit these compounds (Carslaw and Shaw, 2022; Sarwar et al., 2004; Singer et al., 2006; Steinemann, 2016; McDonald et al., 2018). In parallel, this cluster includes C₁₂H₂₇O₄P, corresponding to tributyl phosphate (L4), which does not exhibit a similar seasonal pattern at the rural site. Tributyl phosphate is commonly used as a flame retardant additive, solvent, plasticizer, and antifoaming agent. It has also been reported in aerosol samples collected in indoor environments, such as residential settings and transportation microenvironments (García et al., 2007; He et al., 2018). The co-occurrence of this mixture of compounds suggests a contribution from volatile chemical products (VCPs) and, more broadly, from compounds emitted in indoor environments that undergo transformation and diffusion into the urban atmosphere. In this respect, it noteworthy that recently Zhou et al. (2025) demonstrated that the sources of volatile organic compounds in the atmosphere are rapidly changing showing that indoor environments contribute significantly to the levels of outdoor volatile compounds implying that indoor-to-outdoor emissions significantly contribute to urban VCP levels, rivaling traditional urban sources, e.g., power generation and biomass burning. Supporting the hypothesis of SOA from VCPs, we find in addition 29 compounds attributable to the HOA class of which only 8 were detected at the rural site. The BSOA clusters from both sites share only one HOA compound (C₁₈H₂₉NO₃). However, at SKI, 11 out of 17 HOA compounds were unknown at the urban site. Hence, HOAs contribute significantly to the chemical variability between the two locations.

3.2.3 Agriculture and remaining clusters

HCA revealed additional clusters, which do not exhibit a distinct seasonal pattern. At the rural site we identified a cluster, which display a strong increase during the spring season, with a largest values in May. This is mainly due to a group of molecules belonging to CHNO family and that is responsible for the increase in the overall nTSI of this family shown in Fig. 1. Among these we found C₁₂H₂₅NO, which matches with N,N-dimethyl capramide, used as an additive in pesticide formulations (L4, ECHA, 2024), the aforementioned C₁₂H₂₃NO, attributed to n-octyl pyrrolidone, used in products as solvent for pesticides (L4, ECHA, 2024), C₉H₁₈N₂O, attributed to 1-[[(2R)-oxolan-2-yl]methyl]piperazine (L4), another plant protection agent. By means of standard comparison, we confirmed the presence in SKI samples and the inclusion in this cluster of Metalaxyl (C₁₅H₂₁NO₄, L1, fungicide), Prothioconazole (C₁₄H₁₅Cl₂N₃O, L1, fungicide), Metolachlor (C₁₅H₂₂ClNO₂, L1, herbicide), Pyrimethanil (C₁₂H₁₃N₃, L1, fungicide), Cyprodinil (C₁₄H₁₅N₃, L1, fungicide), Clomazone (C₁₂H₁₄ClNO₂, L1, herbicide), Terbuthylazine (C₉H₁₆ClN₅, L1, herbicide), Dimethenamid (C₁₂H₁₈ClNO₂S, L1, herbicide), and Tebuconazole (C₁₆H₂₂ClN₃O, L1, fungicide). Although some compounds are from verified biogenic (such as pinic and pinonic acids, L1) or probably from other sources (such as C₁₅H₂₄O₄ and C₁₂H₂₃NO, from agricultural residue burning, L4, Zhang et al., 2023a), we labeled this cluster “SKI-Agricultural activities”, pointing out the exclusiveness of many of these features found at SKI. In addition, the cluster fingerprint (Fig. S19) is enriched in S-containing compounds (approximately 14 % of the number of compounds populating SKI-Agricultural activities), of which 34 % are not detected at the urban site. Comparison of the chemical formulas of these S-containing compounds with the results reported by Blair et al. (2017) indicates that approximately 55 % yield a match (L4). Given that Blair et al. (2017) showed that organosulfates can form through the photooxidation of hydrocarbons emitted from diesel and biodiesel in the presence of SO₂, the occurrence of organosulfates, such as ${\mathrm{C}}_{\mathrm{8}}{\mathrm{H}}_{\mathrm{16}}{\mathrm{O}}_{\mathrm{8}}\mathrm{S},{\mathrm{C}}_{\mathrm{16}}{\mathrm{H}}_{\mathrm{26}}{\mathrm{O}}_{\mathrm{3}}\mathrm{S},{\mathrm{C}}_{\mathrm{7}}{\mathrm{H}}_{\mathrm{12}}{\mathrm{O}}_{\mathrm{7}}\mathrm{S}$ (L4), exhibiting temporal trends consistent with those of pesticides and other agricultural activities may suggest a contribution from emissions associated with agricultural machinery (e.g., tractors). However, alternative sources and formation pathways cannot be excluded, highlighting the need to expand available libraries, such as the Aerosolomic database, in order to more confidently link compounds to specific sources. Out of a total of 1398 features that make up the cluster SKI-Agricultural activities, as many as 645 features are completely absent at the Milan site, of which more than 64 % are detected in positive mode. At the MI site, we isolated a cluster populated by CHNO compounds only, mainly measured in positive mode. These features show low nTSI variability along the year with the exception of few events, which show an increase in intensity by orders of magnitude. The fingerprint (Fig. S13) displays that the compounds forming this cluster fall in a very small area of the Van Krevelen diagram ($\mathrm{O}/\mathrm{C}$ < 0.5 and 1.4 < $\mathrm{H}/\mathrm{C}$ < 2.1). Moreover, observing the relation between the chromatographic retention time and the molecular mass, we infer that a subset of these features is likely due to dimers from high-NOx terpene oxidation (Thoma et al., 2025). The median #C (16) value is the highest among the clusters found at MI, while the median OS_c (−1.3) value is the lowest: this is due to the high amount of HOAs, which account for 39 % of the features which belong to this cluster. Very few formulae match with ones found in literature: some of them were attributed to cigarette smoke and cooking by Qi et al. (2019) and to agriculture residue burning by Zhang et al. (2023a). Finally, due to the area span in the Van Krevelen space, and the $\mathrm{O}/\mathrm{N}$ ratio above 3 (3.1 ± 1.4, average ± standard deviation), we labeled this cluster “MI-Aliphatic-OrgNO₃”. The seventh cluster found in MI exhibits a sharp decrease in nTSI in mid-May 2021, stabilizing at nearly constant levels for the remainder of the year. Compounds detected in positive mode outnumber those detected in negative mode. The cluster fingerprint (Fig. S14) shows compounds mainly populating the Van Krevelen space of aliphatics, with a $\mathrm{O}/\mathrm{C}$ ratio close to 0.5 and an $\mathrm{O}/\mathrm{N}$ ratio (for CHNO) of 2.7 ± 0.6, which does not fully restrict the possibilities of N-containing functional groups other than nitrate ester ones. Marker compound matches do not strongly support any hypothesis about the origin of these molecules. We found p-coumaric acid (C₉H₈O₃, L1), which is known to be a product from biomass burning, especially, but not only, from grasses one (Oros et al., 2006). Moreover, C₁₁H₂₂O₅ was found in Li et al. (2021) (L4) during oxidation of n-dodecane experiments and C₇H₁₃NO₂ was detected as m-xylene SOA (L4, Li et al., 2018). Due to the sharp decrease of the nTSI we can speculate that the driving source of this cluster is an anthropogenic activity but no clear hint is found. Thus, we labeled this cluster “Other anthr-OA”.

3.3 Organic aerosol molecular variability between Milan and Schivenoglia

The use of the HCA technique allowed us to characterize cluster of compounds with similar temporal pattern. Between the two sampling sites, similarities and differences in molecular composition of the interpreted clusters are discussed in Sect. 3.2. Their contribution to the annual nTSI is also diverse: winter-peaking clusters in Milan, linked to traffic, biomass burning, and a mix of both, account for about 16.1 %, 10.9 % and 8.1 % to annual nTSI, respectively. At Schivenoglia, mixed traffic and biomass burning sources account for 14.9 %, the cluster linked to biomass burning with ammonium nitrate formation account for 20.2 %, and the third cluster containing aliphatic amines possibly from agricultural reside burning, accounts for 2.6 % to the annual nTSI. The compounds attributed to biogenic-related sources, peaking during summertime, have a higher impact at MI (31.4 %) rather than at SKI (26.4 %): in addition, a fraction of MI-Summer B/ASOA cluster (11.0 %) might also be due to biogenic sources or stem from VCP (e.g. limonene) oxidation. At the urban site, anthropogenic sources have an additional 14.1 % contribution to the annual nTSI due to MI-other anthr-OA cluster, and 2.1 % due to MI-Aliphatic OrgNO₃ one, even though no clear attribution is performed. At the rural site, anthropogenic activities linked to agricultural practices are clearly observed and quantified to 27.8 %. Notably, some clusters of MI and SKI show a similar molecular composition. For example, 61 % of compounds clustered in MI-BSOA are in common with SKI-BSOA, and in turn SKI-BSOA shares 56 % of its compounds with MI-BSOA and 12 % with MI-Summer B/ASOA. Both biogenic clusters show a number of compounds that are not classified at the other site: at MI 350 compounds are not detected (or detected in less than 10 % of the analyzed samples) in SKI, whereas a higher amount (436) are observed at the rural site but not at the urban one. This provides a first insight of the molecular variability between the two environments. For SKI, the cluster that contains the largest number of these site-specific compounds is SKI-Agricultural activities (645, i.e. 46 % of compounds falling into it). On the contrary, every cluster found in MI contains a comparable amount of site-specific compounds. The chord chart (adapted from Liu, 2025) in Fig. 4, graphically assists the assessment of the site differences, where the width of the flow-lines is proportional to the number of compounds that are common between the two connected clusters. The chart suggests that combustion sources have a different profile at the rural compared to the urban site. MI-BBOA, MI-TrOA and MI-Other combustion share 45 % (n=412), 19 % (n=175), and 33 % (n=306) compounds only with SKI-BBOA/TrOA, and 11 %, 30 % and 15 % with SKI-Nitro Aromatics, respectively. This leads to a percentage of site-specific combustion-related compounds from this MI clusters of 36 %, 40 %, and 30 % respectively.

Figure 4 Chord chart illustrating, for each cluster, the number of compounds shared between the two sites, as well as the number of cluster-specific compounds unassigned at the other site, which are considered site-specific.

We found a total of 2513 and 1525 compounds that are site-specific for MI and SKI, respectively, confirming the higher molecular-diversity at the urban site. The site-specific fingerprints, shown in Figs. S20 and S21, highlight that the urban site is richer in nonpolar compounds with high molecular mass, and in CHO-containing aromatics. Among aliphatics, N-containing compounds are the most prominent compounds present at the urban site. Figure 5 illustrates the contribution to the monthly nTSI of each cluster at both sites. As expected, during wintertime the nTSI explained by the site-specific compounds (in darker colors) is lower than in the other seasons. The reason is explainable by the lack of both weak oxidant conditions, leading to a longer lifetime of organic molecules, and meteorological events able to abate PM concentrations: as a result, aerosols disperse, mix and age in the lower layers of the troposphere in a vast area of the Po basin. Limiting the discussion on the site-specific features only, during winter in both sampling sites the anthropogenic clusters explain between 10 % and 15 % of the nTSI, whereas site-specific BSOA compounds are almost negligible. On the opposite, from May to September this cluster provides the highest molecular variability at MI in respect of SKI. The rural site, in fact, exhibits site-specific compound profiles predominantly originating from agricultural activities throughout the year, with only modest contributions from the SKI-Nitro-Aromatics cluster in winter and from source SKI-BSOA in summer. Between May and September, both sampling sites exhibit a significant molecular variability, which rises from 15 % to 25 %. It is worth noting that, at both sites, the lowest concentrations of OC were measured in May (2.0 ± 0.6 and 1.6 ± 0.3 µg m⁻³, Fig. 1). At the same time, the highest plant protection product intensities at the SKI site were detected in May, suggesting that, despite theoretically indicating a reduced health impact that month, the aerosols could have a non-negligible toxicological effect during that period at SKI. Simultaneously, the MI-Summer B/ASOA cluster, predominantly composed of compounds detected in positive mode, with a significant fraction unknown at the SKI site, shows its highest contribution to nTSI (7 % accounting for site-specifice compounds only) once again in May. These results highlight the need to assess the health effects of aerosols not only based on the mass quantification metrics but also, and perhaps most importantly, based on chemical composition and evaluate health effects of emissions from agricultural practices.

Figure 5 Contribution to the normalized total signal intensity (nTSI) of the classified compounds divided by cluster attribution. Darker colors represent the fraction of site-specific features of each cluster.

Figure 6 Cluster attribution of the contributions to the nTSI from residuals, limited to pairs of samples that show similar OC concentrations and OC contribution to PM_2.5. Darker colors indicate the fractional contribution from compounds that are site-specific to each sample.

Finally, we select sample pairs with similar OC concentrations (within a tolerance of 2.0 µg m⁻³) and contribution to PM_2.5 concentration (within a tolerance of 5 %) collected on the same day at both sites (hereinafter referred as “similar samples”). Considering the full year, 116 sample pairs meet this criterion, although they are unevenly distributed over time (Fig. S22). In particular, the number of pairs is the lowest in June (#5), September (#1), and December (#1). When restricting the analysis to the subset of filters used in this study, the number of valid pairs decreases to 54, and no pairs meet the above criterion in June and December. Consequently, these months are not represented in the figures showing similar sample pairs (Figs. 6 and 8). These conditions lead us to select 54 sample pairs across all months but June and December 2021. The differences of the nSI highlight the dissimilarities between sample pairs of similar OC concentrations. Once again, we take advantage from the HCA performed (Sect. 3.2) and classify the enrichment in compounds based on their source attribution. The bar chart in Fig. 6 shows the compound membership clusters with residuals, averaged on a monthly scale, favorable to MI (upper semi-axis) or SKI (lower semi-axis) for the selected samples as contribution to nTSI (see Sect. S1.2, “Residuals contribution calculation”). Both semi-axes are plotted as absolute values to display only the magnitude of the residuals. At MI site, the molecular variability is mostly due to compounds detected at both sites (lighter colors) for the majority of the sources with a seasonal pattern. The residuals of common compounds contribute about 38 ± 7 % to nTSI (reported as an annual average), whereas the residuals of site-specific compounds contribute 18 ± 5 %. At SKI residual contributions of common compounds account to 26 ± 6 %, slightly higher than for site-specific ones (21 ± 11 %). In particular, Fig. 6 shows that the agricultural activities drive both the site-specific features and the high residual contributions. Moreover, the highest differences are observed in May, when the OC concentrations are the lowest. This reinforces the assertion that the sources at the two sites have distinct emission profiles and formation pathways.

4 Atmospheric implications on light-absorbing aerosols

In the initial dataset, we found 4370 compounds with $\mathrm{DBE}/\mathrm{\#}\mathrm{C}$ ranging between 0.5 and 0.9 (Fig. 2), of which 1551 and 1247 were classified for MI and SKI, respectively, using HCA. Among these, we detect C₈H₇NO₅ attributed to 3-methoxy-4-nitrobenzoic acid (L4, Li et al., 2020a), C₇H₆O₂, attributed to benzoic acid or p-hydroxybenzalehyde (L4, Worton et al., 2011), p-coumaric acid, nitrophenol, nitrocatechol, acetosyringone (L4), phtalic acid, C₁₀H₁₀O₃ assigned to coniferyl aldehyde (L4, Oros et al., 2006 and Huang et al., 2022), C₆H₆N₂O₃ attributed to amino-nitrophenol (L4, literature within Huang et al., 2025), C₆H₃N₃O₇ attributed to picric acid (L4, Lin et al., 2015), C₈H₆N₂O₂, likely nitroindole (L4). Potential BrC compounds are listed as Additional data (Table A2) together with their cluster attributions and the average percentage contribution to nTSI. The fingerprints of the classified BrC features at MI and SKI are shown in Figs. S23 and S24. The contribution of classified BrC-related compound to nTSI accounts for about 6 % at both sites during summer and increases up to 38 % during winter, with an annual average of 14 % at MI and 13 % at SKI, respectively. Most of BrC-related compounds fall into combustion clusters (Fig. S25): at the urban site, we observe that all BrC-combustive clusters exhibit a stronger contribution during wintertime, with the highest one due to the MI-Other combustion cluster, which is mainly populated by aliphatic-CHNO and CHN(+) compounds. MI-TrOA cluster seems to have less impact compared to MI-BBOA, which still shows the greatest contribution for combustion sources during the summer months. At the rural site, however, the largest contribution comes from the “SKI-NitroAromatics” cluster, which is also linked to combustion but, unlike at MI site, appears to have an increasing impact on nTSI until early spring. Here too, SKI-BBOA/TrOA contributes to the presence of potentially BrC compounds throughout the year, though with a smaller contribution relative to nTSI when compared to the role of the respective clusters found at MI. Considering the 54 sample pairs with very similar OC concentrations and contributions to PM_2.5 mass (Sect. 3.2.3), Fig. 7 present nested annular plots illustrating the residual contributions of all (outer circle) and BrC-related compounds only (inner circle), averaged over the entire year. The plots confirm, as expected, that combustion sources are the dominant contributors to the annual variability of BrC, accounting for nearly 75 % of the residual contributions at both sites. However, the average contribution of site-specific BrC-related compounds is notably lower but still relevant: they constitute 23 % of the residuals at the urban site, compared to 16 % at the rural site. During the year the residual contribution of the clusters to nTSI does not vary substantially (between 2 % and 9 % for all classified BrC-related compounds, as shown in Fig. 8, secondary vertical axis) at both sites, even with a seasonality. However, because differences are observed both in abundance (residuals are not null) and in the molecular species present, a different interaction with light can be hypothesized, even under similar OC concentration conditions.

Figure 7 Residuals calculated on samples with similar OC concentrations and contributions to PM_2.5 for MI (upper plot) and SKI (lower plot). The outer circle illustrates the distribution of annual average residuals across all compounds attributed to sources identified at the reference site. In contrast, the inner circle focuses on the same distribution but is limited to compounds potentially classified as LAA. Darker shades highlight the contribution of the site-specific compounds.

Figure 8 BrC residuals at MI (top) and SKI (bottom) for samples with similar OC concentrations and contribution to PM_2.5. Darker colors highlight site-specific BrC features. The secondary vertical axes show the monthly averaged contribution of the residual BrC-related compounds to the nTSI of the selected samples.

In this regards, the AAE of eight sample pairs collected at MI and SKI sites provides an insight of this (Fig. S26). For the majority of samples, the AAEs estimated for MI samples were higher than those collected at SKI: the median absolute differences were 1.40, 0.23, and 1.22 in the 230–250, 250–400, and 400–650 nm ranges, respectively, with absolute maximum up to 4.90. Although more measurements are needed, these results confirm that, even with a low OC concentration variability, the molecular differences between two sites lead to a different light-absorption behavior.

5 Conclusions

This study explores the chemical variability of organic aerosol in a European hotspot for air quality, such as the Po Valley. Historical data indicate that PM_2.5 and OC concentrations are often similar across sampling sites, suggesting little variability in aerosol composition. To test this, PM_2.5 samples collected in 2021 at a rural agricultural and an urban background site were analyzed. Using a hierarchical clustering approach on UHPLC-HESI-UHRMS data, compounds were grouped by known source markers and temporal patterns. At the urban site (Milan), seven clusters emerge: three peak in winter (linked to traffic combustion, biomass burning, and a mix of both), two in summer (one due to biogenic secondary organic compounds, the other made of biogenic and anthropogenic). The remaining two clusters without a clear seasonal trend likely represent unidentified anthropogenic activities. At the agricultural site (Schivenoglia), five clusters are identified, three of which peak in winter. One correlates with soluble potassium, levoglucosan, benzene, and NO₂, suggesting combustion sources; another links biomass burning with ammonium nitrate formation, and the third contains aliphatic amines possibly from agricultural residue burning. A biogenic SOA cluster peaks as expected during summer, and one cluster clearly emerges in late spring tied to agricultural activities, as its main compounds are plant protection products. Focusing only on site-specific compounds of each site, the molecular richness of the urban site amounts to about 38 % of the total classified features, contributing between 11 % and 23 % to the nTSI on monthly basis. The agricultural site is characterized by 27 % site-specific features, contributing between 7 % and 23 % to the nTSI, with a peak of up to 55 % in May due to molecules associated with plant protection products. At both sites, the lowest contribution from molecular variability is observed in the winter months, suggesting that poor ventilation of the basin due to meteorological conditions results in a more homogeneous Po Valley during this season. The results indicate that, while the two sites may appear quantitatively similar, they exhibit significant qualitative differences. Although limitations in the analytical approach may influence the findings, the observed variability suggests that it could distinctly affect both the toxicity of organic aerosols and their interactions with radiation. In this regards, we demonstrate that PM_2.5 extracts have significant spectral behavior differences. This work encourages the aerosol research community to deepen the molecular-level understanding of organic aerosols and to identify relationships with larger-scale effects, thereby guiding policy makers toward the adoption of more effective measures. Finally, the use of high-resolution and high-accuracy mass spectrometry enables the detection and identification of specific compounds based on targeted hypotheses and data contextualization. This approach allowed us to identify the contributions of certain sources, such as agricultural activities, which are often overlooked in traditional source apportionment methods. Although our analysis cannot be quantitative, it provides a more refined understanding of the contributions of various anthropogenic sources.