Interlaboratory comparison of humic substances compositional space as measured by Fourier transform ion cyclotron resonance mass spectrometry (IUPAC Technical Report)

2.1 Intercomparison set of humic materials and FTICR MS facilities

Six HS samples were included in the experimental set: one sample of humic acids (HA) from the Suwannee River (SRHA, the IHSS standard), one sample of humic acids from the sod-podzolic soil collected near Moscow (SHA-PwM, the state standard sample, Russia), one sample of humic acids from the Mollisol collected near Kursk (SHA-CtK, the state standard sample, Russia), one sample of fulvic acids from the sod-podzolic soil collected near Moscow (SFA-PwM), one sample of fulvic acids from the Mollisol collected near Kursk (SFA-CtK), and one sample of coal (leonardite) humic acids (CHA) isolated from commercial humate. The samples were provided by the laboratory of Dr. Irina Perminova (Lomonosov MSU, Moscow, Russia). The selected samples were delivered to five research groups who had access to FTICR mass spectrometers. These five groups are listed in Table S1.

2.2 The guidelines for intercomparison study

2.3 Data collection and formulae assignment

The obtained FTICR MS results were collected in the form of calibrated non-filtered mass lists with signal to noise ratio exceeding four and anonymized by assigning codes to each laboratory in randomized order (a, b, c, d, e) prior to data treatment. The raw mass lists are provided in the Supporting Information. All data were treated for formula assignment using house-made Transhumus software based on the open-source R environment designed by A. Grigoryev. The parameters used for formula assignments were set with sensible chemical constraints according to the literature: O/C ratio below or equal to 1, H/C between 0.3 and 2.2, element counts (C atoms up to 120, H atoms up to 200, oxygen atoms between 1 and 60, N atoms up to 2, and sulfur atoms up to 1), and a mass accuracy window below 0.5 ppm (Da/MDa) [19]. The data analysis was performed by the Lomonosov MSU research group shown in Table S1.

2.4 Statistical treatment

For direct comparison of the mass spectra acquired by the different research groups for the same HS sample, we calculated a Jaccard similarity score (J score) based on the presence/absence of the formulae in the pairs of the corresponding mass spectra [20]. The J score is formally defined as the size of the intersection divided by the size of the union of the sample sets. The mathematical representation of the J score is written as:

where the numerator is the number of common formulae in the mass spectra measured for the same HS sample (area of overlap), and the denominator is the total number of formulae found in the different spectra for the same HS sample (area of union).

For correlation analysis, the FTICR MS assignments were plotted in a van Krevelen diagram, where the x and y axes represent a ratio of the numbers (n) of O to C atoms (nO/nC) and a ratio of the numbers of H to C atoms (nH/nC), respectively, as described by Kim et al.[21]. The obtained diagram was quantitatively treated using a cell-partitioning algorithm by binning the diagram area into 20 cells, as described by Perminova [22]. The partitioning is shown in Fig. 1, using the SRHA sample as an example. The cell-based distribution of experimental points was calculated by quantifying the intensity-weighted population density of each distinguished squared area in a van Krevelen plot cell (D_i), as expressed by Equation (2) below:

where D_k is the intensity weighted population density of the cell k and N is the total number of points in the van Krevelen diagram. In the numerator, N_k is the number of points belonging to the particular cell k and I_i is the relative intensity of a peak from the mass-spectrum represented as a point i within the cell k. while in the denominator I_j is the summarized relative intensity of all peaks of the mass-spectrum corresponding to the points populating in the diagram.

Fig. 1:

An example of binning a van Krevelen diagram into 20 cells [22], which was used for data treatment in this study. An example is shown for the SRHA sample.

The obtained values of population density could be used as descriptors of the chemical space of HS and NOM, per se, as 20 values – this approach was used in this study. They can be also combined based on the classes of chemical compounds, using the mapping of the van Krevelen diagram as described in [22]. Another approach to splitting the van Krevelen diagram is based on a use of the aromaticity index (AI) [23]. Atomic constraints for the AI-based partitioning were proposed by Kellerman et al. [12], as is shown in Table 1.

3.1 Intercomparison of the FTICR MS data for the SRHA sample

We started by examining the mass spectra obtained for the IHSS standard sample, SRHA, by five participating laboratories (Figs. 2 and 3). Mass distributions turned out to be substantially different, in spite of the same protocol for sample preparation and the same given conditions for instrumental setup being used by the participants. The spectra of the laboratories a, b, and c contained very intense peaks of impurities. For example, mass-spectrum obtained in laboratory a was dominated by a low-mass even peak, whereas SRHA is composed mostly of odd peaks. This could be indicative of plastic and fatty acid impurities. At the same time, the mass spectra obtained in laboratories d and e did not have peaks of impurities. Nevertheless, due to high sensitivity and the dynamic range of FTICR mass-spectrometers, the number of peaks with a signal to noise (S/N) ratio larger than four was about 20 000 in all cases, except for laboratory c, which showed the lowest number of peaks. The difference in mass spectra was also reflected in peak distribution at nominal mass, as is shown in Fig. 2. For all five spectra, we observed more than 10 monoisotopic peaks at a nominal mass. The maximum intensity for laboratory a was observed for the ion with m/z equal to 301.05651, which corresponds to the unsaturated compound with molecular formula C₁₂H₁₃O₉ (H/C equals to 1.16). At the same time, for laboratories a, b, c, and d, the maximum intensity was attributed to the ion with m/z equal to 301.03538, belonging to an aromatic compound with molecular formula C₁₅H₉O₇ (H/C equals to 0.66). Abundant saturated ions with high mass defect were observed only in the case of laboratory c, while they were absent in the other spectra, as is shown in Fig. 3.

Fig. 2:

FTICR mass spectra of SRHA obtained by five participating laboratories. Insets show the magnified mass-spectra fragments for a nominal value of m/z of 301.

Fig. 3:

Van Krevelen diagrams plotted for CHNOS assignments made for the SRHA sample from FTICR MS data obtained by the five intercomparison laboratories: a, b, c, d, and e.

We assigned all molecular formulae for monoisotopic peaks, excluding heavy isotopologues, such as ¹³C and ³⁴S. The general characteristics of the data set obtained for SRHA are given in Table 2. In the obtained range of data, laboratory a determined the highest number of assigned molecular compositions for SRHA (6403), while laboratory d determined the lowest (1748). Despite substantial variation in the number of assigned formulae, all five spectra were characterized by the dominant contribution of CHO compositions. The intensity of all other stoichiometries did not exceed 12 %. Very substantial variations were observed in the number of CHOS formulae. Where laboratories a and c determined 499 and 213 compositions, respectively, laboratories b, d, and e found significantly fewer of these compositions, lower by a factor of 10. The abundant CHOS compounds could also be originated from impurities, e.g., surfactants [24].

Despite the differences in the number and abundance of the assigned formulae, the calculated number-averaged parameters, such as M_n, AI, nO/nC, and nH/nC, revealed similarities between the spectra examined. Four out of five spectra (except for those of laboratory d) were characterized by number-averaged AI (AI_n) of about 0.5 and M_n between 450 and 500 Da. Values of AI_n equal to 0.3 and M_n equal to 379 Da were observed for laboratory d, which indicates the higher contribution of low molecular mass and relatively saturated compounds, as compared to the other laboratories. The number-averaged values of O/C ratio also indicate higher contributions of reduced compounds in laboratories b and d, with nO/nC below 0.5, compared to laboratories a, c, and e. This could be explained by the different parameters of ion accumulation in ion source and collision cell adjusted in experiments, which are known to affect the proportions of oxidized and reduced components [25].

We plotted the obtained data into a van Krevelen diagram [21] to visualize the differences in the molecular compositions determined by five different laboratories (Fig. 3).

The van Krevelen diagrams obtained for the same sample were characterized by similar general features of molecular distribution, but with substantial differences in the exact location of molecular species. The SRHA spectra measured by the b and d laboratories were characterized with the highest abundance of the reduced compounds. These clear outliers affected the number-averaged values given in Table 2. The data obtained by the a, c, and e laboratories were characterized by the dominant contribution of oxidized aromatic species with an O/C atomic ratio between 0.5 and 0.7, which is inherent for NOM from temperate and tropical regions [26], [27]. The data sets of SRHA-b and SRHA-d did not show any distinct features. At the same time, the SRHA-d was fully depleted of oxidized compounds and its whole molecular profile was shifted towards saturated compounds.

The common and unique stoichiometries were calculated and plotted in van Krevelen diagrams, as shown in Fig. 4a and b, respectively. In total, there were 9193 formulae determined for the SRHA sample by five spectrometers. Of these, 1260 were “common” (determined by all five spectrometers), 4272 were “shared” (detected by two or more spectrometers), and 3661 were “unique” (detected by only one spectrometer). Fig. 4a shows that the stoichiometries with the most robust values, which were determined in the SRHA sample by all instruments, were located in the center of the van Krevelen diagram, close to the O/C value of 0.5 and H/C value of 1. Another observation was that the least robust stoichiometries (Fig. 5b) were positioned in the bottom left corner of the van Krevelen diagram, which is usually prescribed to “black carbon” [28]. As a rule, these components possess the lowest abundance in FTICR mass spectra; they are poorly resolved due to low S/N ratios. This produces a large uncertainty in formulae assignment due to variability in the sensitivity and dynamic range of different instruments.

Fig. 4:

The van Krevelen diagrams for the common (A) and unique (B) stoichiometries determined in the SRHA sample by five different FTICR MS instruments.

Fig. 5:

Compound class distributions in the mass spectra of the SRHA sample as measured by five participating laboratories (x-axis).

The example of direct comparison of individual stoichiometries obtained for the same sample (SRHA) by five spectrometers is a good visualization of very small overlap (14 %) in the discrete stoichiometries of HS samples determined by different FTICR MS instruments. Given much better similarity in the intensity-weighted parameters, shown in Table 2 and similar molecular profiles in Fig. 4, we have binned the area of the van Krevelen diagram into nine regions, assigned to different compound classes, as is described in Table 1, and calculated the occupation density of each area. The obtained data are shown in Fig. 5

All spectra showed the prevalence of unsaturated compounds in the molecular composition of SRHA, which provided for 60 % of the total intensity of formula assignments. Still, the spectra obtained by the b, d, and e laboratories were characterized by the prevalence of low-oxidized unsaturated compounds, while the spectra acquired by the a and c laboratories were characterized by the major contribution of highly-oxidized unsaturated compounds and a similar abundance of low- and highly-oxidized species, respectively. The d laboratory determined an extremely high content of saturated compounds, which can be considered the outlier. For the other four laboratories, the ratio of the summarized relative abundance of unsaturated compounds to aromatics was two, whereas for the d laboratory it was equal to 3. Thus, all participating laboratories revealed the major role of substituted aromatic structures indicating the ligninic nature of SRHA. The ratio of oxidized to reduced species varied among the laboratories, which can be explained by different equipment adjustments for analysis. The obtained results highlighted the importance of the development of statistical tools for data treatment, which could explore the unique and common features of HS samples analyzed by different FTICR MS instruments.

3.2 General characteristics of the full data set on intercomparison of FTICR mass spectra

The collected data set of FTICR MS data was composed of 30 spectra: five spectra for six HS samples measured by five participating laboratories. The raw mass-spectrometric data were compared using two signal properties: m/z values and relative abundance, expressed as a logarithm of signal intensity for better visualization. We highlighted the “shared” peaks in blue (“shared” means the peaks observed by two or more laboratories with an uncertainty of less than 0.2 ppm) and the “unique” peaks (i.e., observed by a single laboratory) were highlighted in orange, as is shown in Fig. 6. It would be reasonable to assume that the signals with the highest intensity could be observed in all the spectra of one sample, and the low-intensity ones would be those unique to different spectra. In that case, it would be possible for each mass range to find a threshold separating the common signals from the unique ones. It would also be logical to assume that such a threshold would depend on the mass: it could increase or decrease. Under these circumstances, one could find such a straight line that would separate a set of the common (intense) signals from the set of unique (low-intensity) ones. However, in this study it was shown that such separability was not observed within the data. Hence, our initial assumption was wrong.

Fig. 6:

Thirty FTICR mass spectra (six samples and five laboratories) collected for the intercomparison study. Relative abundance is represented in logarithmic scale. Shared (present in at least two spectra) and unique peaks are highlighted in blue and orange colors, respectively.

This might indicate that there is no simple rule for discrimination between the “common” and the “residual” formulae. Still, abundance (or peak intensity) had substantial discriminating power with regard to the common and residual peaks. For the four laboratories a, b, c, and d, the shared peaks were much more abundant than the unique ones in four samples out of six (this was not the case for those determined by the e laboratory). The highest discrepancy between the measurements was observed for the SHA-CtK and CHA samples, which were characterized with the largest number of abundant unique peaks. It is important to note that these very samples were characterized with the highest contribution of condensed aromatic compounds [29]. Taking into account the low ionization efficiency of these molecular constituents of humic ensemble, we suggest that the number of these species might be critical when comparing results from equipment with different analytical characteristics, as in the case of this work. The particular problems can be seen in ionizing the samples enriched with condensed components for laboratory d: the mass spectra in this case had a few common peaks with those observed by the other laboratories.

The high number of unique peaks detected by different laboratories might be explained by three reasons: discrepancy in a mass measurement, instrument-specific ion accumulation parameters, and different dynamic ranges and sensitivities of the equipment used. Recommendations for mass-spectra acquisition and calibration described in Section 2.2 of this manuscript were aimed at minimizing the first parameter. However, the dynamic range and sensitivity are specific features of the instruments and depend mainly on the ion optics, ICR cell geometry, and the homogeneity of the magnetic field [30]. In addition, the lifetime of ions depends drastically on the ion accumulation parameters [24]. The latter were not standardized in this study due to the different architecture of the instruments used by the participating laboratories. In Fig. 6, it can be seen that even after peaks filtration by S/N equal to 4, the laboratory a possessed the widest dynamic range at 10⁸ arbitrary units of peak intensities, while the narrowest dynamic range, at 10⁵, was observed for laboratory e.

The next step was a comparison of the CHNOS formulae assignments, which were made for the obtained FTICR MS data sets. The data are shown for the total number of the assigned formulae and for the number of unique formulae in Fig. 7a and b, respectively. From Fig. 7, it can be concluded that the initial data acquired by laboratories differed substantially in the total number of assigned formulae for the intercomparison HS sample set. The lowest number of formulae was assigned to the soil HA from Mollisol (SHA-CtK). The number varied from as low as 1743 (laboratory b) to as high as 3380 (laboratory e).

Fig. 7:

Number of the formulae that could be unambiguously assigned to six intercomparison samples based on the FTICR MS data of the five participating laboratories. A) the total number of assigned formulae; B) the number of unique formulae. The color bar shows the number of formulae, from the darker blue color, representing fewer formulae, to the lighter colors, representing more formulae.

The total number of assigned formulae was followed by three other HA samples: SRHA, CHA, and SHA-PwM. A greater number of assigned formulae were observed for the more oxidized and lower molecular weight soil fulvic acids: in particular, sod-podzolic soil FA (SFA-PwM). With respect to the laboratories, the highest values were provided by laboratory a (the values ranged from 3000 to more than 10  000), followed by laboratory e (the values ranged from 3000 to 9000). Laboratories b and d provided lower values, from 1700 up to 7000, while laboratory c provided the lowest numbers, from 1700 to 3600. It was also of interest to examine the number of unique identifications made by the different laboratories (Fig. 8b). In the cases of laboratories a, b, and c, the number of unique identifications was higher for FA than HA, while an inverse pattern was seen in laboratories d and e. Laboratory a provided the highest number of unique formulae for all samples, while laboratories b and c yielded the lowest number of unique identifications. The dramatic discrepancy between the number of detected peaks and the assigned formulae clearly illustrates the challenge of direct comparison between HS molecular compositions determined using different instruments.

Fig. 8:

The occupation density van Krevelen diagrams plotted from the CHONS formulae assignments for the six samples used in this study (from top to bottom): CHA, SHA-CtK, SHA-Pw, SFA-CtK, SFA-Pw, and SRHA, based on the results obtained by five intercomparison laboratories (from left to right): a, b, c, d, e.

To compare molecular profiles determined by FTICR MS for six HS samples under study, all C_cH_hN_nO_oS_s assignments were used for the calculation of H/C and O/C ratios and plotted in van Krevelen diagrams as population density maps versus the common representation using discrete stoichiometries (Fig. 8). This type of data presentation allows for robust fingerprinting of integral common features in the formulae distribution over the 2D space [31].

For all HS samples studied, the substantial similarity of occupation density distribution was observed for the same sample measured by different laboratories. In the case of CHA and SHA-CtK, the aromatic and condensed compounds possessed the highest abundance. At the same time, the SHA-Pw sample was characterized by the significant contribution of saturated compounds. Both SFA samples and the SRHA sample were characterized by the high contribution of oxidized aromatic compounds. These trends are in sync with the data on the elemental and structural group composition analysis of the samples used in this study [32], [33], [34]. Despite the clear outliers obtained for CHA and SHA-Ctk by laboratory d, the instruments were capable of discrimination between the HS fractions used in this study.

3.3 Similarity analysis of the FTICR MS data obtained in this study

Statistical analysis was used for assessing the similarities between the molecular compositions determined for the six HS samples used in this study by the five participating FTICR MS laboratories. Similarity heatmaps were calculated for all pairs of analyzed samples based on Euclidian distances between vectors in multidimensional space composed of the abundances of the assigned molecular formulae. The logic behind this approach is that the J score can be calculated for any two spectra (see Equation (1)). Then, the matrix with a size of N × N can be composed out of N spectra, where the positions (i, j) will be taken by the J score values varying from 0 to 1. The value 0 indicates an absence of common signals, whereas the value one indicates full coincidence of the spectra (up to differences in intensity). It is obvious that such a matrix will have one values on the diagonal, and that it will be symmetric. To graphically display such a matrix, the values could be color-coded and designated in the color bar. This type of matrix is called a “heatmap”. The type of heatmap depends on the order of spectra arrangement. If the spectra are grouped according to certain criterion, for example, according to the sample type (six samples measured by five laboratories), then the heatmap will show a 5 × 5 submatrix that characterizes the similarity of each sample across all laboratories. The heatmaps constructed for the data obtained in this study are shown in Figs. 9 and 10.

Fig. 9:

Similarity heatmaps for all pairs of the analyzed samples based on Euclidian distances between vectors composed of molecular formulae abundances; A) the heatmap for samples sorted by the source; B) the heatmap for samples sorted by the laboratory.

Fig. 10:

The values of Jaccard similarity scores calculated for the presence/absence of the peaks in the initial FTICR mass spectra for the six HS samples measured by five laboratories. The color bar refers to a value of the J scores.

The grouping of the samples was made by source and by laboratory (Fig. 9a and b, respectively). We observed poor similarity between the data obtained by different laboratories for all HA samples (Fig. 9a). Only in case of SHA-Pw, the certain degree of similarity was observed between the laboratories b and c. At the same time, both SFA samples and the SRHA sample were characterized by significant similarity between the labs. This result could be related to the much better ionization efficiency of the polar components of FA compared to more aromatic and less oxidized HA, yielding much more comparable results.

The second tool that we used for similarity analysis was the comparison of the presence/absence of the formulae in the data set using the Jaccard similarity score calculation. The results are shown as heatmaps in Fig. 10. As can be deduced from Fig. 10, the analysis for the presence of the same peaks in the spectra measured by the different laboratories yielded poor results, with very low J score values, i.e., not exceeding 0.56 (not more than 56 % of the pairwise shared peaks). The obtained values can be interpreted as meaning that, for each peak in the spectrum, there is no confidence of its presence in the spectrum obtained by a different laboratory.

In search of ways to improve the similarity between the data obtained for the same sample by different laboratories, we have generated descriptors for the density distribution of the assigned formulae over the field of the van Krevelen diagram. For this purpose, we binned a van Krevelen diagram into 20 cells within the range of nH/nC values, from 0.2 to 2.2, and nO/nC values, from 0 to 1, as is shown in Fig. 1, followed by a calculation of the intensity-weighted contribution of each cell (D_k) using Equation (2), as proposed by Perminova [22]. Then, we conducted a correlation analysis using the obtained descriptors and calculated the squared value of correlation coefficients (R²) for each pair of mass spectra (i.e., laboratories) for each sample. The results of this correlation analysis are shown in Fig. 11.

Fig. 11:

The results of pairwise correlations for the pairs of laboratories for each sample (A) and for the pairs of samples for each laboratory (B) undertaken for the population density of a van Krevelen diagram binned into 20 cells and described using population density descriptors (D_k,Equation (2)), as proposed by Perminova [22]. The color bar refers to R² values.

This analysis has shown the presence of two outliers in the data set, nominally, the CHA and SHA-CtK spectra measured by laboratory d. These were characterized with very low R² values, seen as dark crosses in the corresponding matrices in Fig. 11a. This is consistent with a visual inspection of the van Krevelen diagrams of these samples, which differ substantially between the laboratories. Except for these clear outliers, all other correlations were characterized with high R² values (close to 1), which shows a good applicability of this approach to the similarity analysis. It was also of interest to explore the correlations between the samples or to run the correlation analysis for each pair of samples from each laboratory. The corresponding data are shown in Fig. 11b. Except for laboratory d, we observed very high R² values for the samples with similar origin, for example for both SFA samples and the SRHA sample, two SHA samples, and the coal HA sample, seen as the bright green and yellow squares in Fig. 11b. The correlation analysis indicates that the different FTICR MS instruments acquired similar results for the density distribution of the assigned formula over the van Krevelen diagram for the HS samples from different sources used in this study.

While conducting the above data analysis, we have established that the presence of the zero values corresponding to the non-populated cells in the correlation data set increased the correlation relationship between the sets on the account of “0” data. These “0” data can be clearly seen in Fig. 1, representing occupation density D_k (Equation (2)) of the least populated cells, e.g., the cells from 1 to 5, or from 15 to 20. This is why we ran the corresponding analysis by excluding all 0 values (i.e., non-populated cells) from the data set. A good countermeasure for the appearance of 0 values is to enlarge the areas of the van Krevelen diagram. This can be done by ascribing a sum of the D_k values to the certain compound class using mapping as described by Perminova [22]. Another approach is to bin the van Krevelen diagram with regard to the AI value, as described by Koch et al. [13]. We conducted this alternative binning of the van Krevelen diagram into the nine areas, as described in Table 1, and used the obtained descriptors for the correlation analysis. The results are shown in Fig. 12.

Fig. 12:

The results of pairwise correlations for the pairs of laboratories for each sample (A) and for the pairs of samples for each laboratory (B) undertaken for the population density of a van Krevelen diagram binned into nine areas in accordance with the AI value (Table 1) as proposed by Kellermann et al. [12]. The color bar refers to the R² values.

It was found that a use of the coarser binning of the van Krevelen diagram into nine areas gave slightly worse results than the finer splitting into 20 cells, both for the pairwise correlation (each pair of laboratories for each sample) and when pairing by laboratory (each pair of samples for each laboratory). However, the results were in general consistent and supportive of a high similarity of population density of the compositional space of the six HS samples used in this study measured by five FTICR MS spectrometers.