Real-Time Volatile Metabolomics Analysis of Dendritic Cells

Dendritic cells (DCs) actively sample and present antigen to cells of the adaptive immune system and are thus vital for successful immune control and memory formation. Immune cell metabolism and function are tightly interlinked, and a better understanding of this interaction offers potential to develop immunomodulatory strategies. However, current approaches for assessing the immune cell metabolome are often limited by end-point measurements, may involve laborious sample preparation, and may lack unbiased, temporal resolution of the metabolome. In this study, we present a novel setup coupled to a secondary electrospray ionization-high resolution mass spectrometric (SESI-HRMS) platform allowing headspace analysis of immature and activated DCs in real-time with minimal sample preparation and intervention, with high technical reproducibility and potential for automation. Distinct metabolic signatures of DCs treated with different supernatants (SNs) of bacterial cultures were detected during real-time analyses over 6 h compared to their respective controls (SN only). Furthermore, the technique allowed for the detection of 13C-incorporation into volatile metabolites, opening the possibility for real-time tracing of metabolic pathways in DCs. Moreover, differences in the metabolic profile of naïve and activated DCs were discovered, and pathway-enrichment analysis revealed three significantly altered pathways, including the TCA cycle, α-linolenic acid metabolism, and valine, leucine, and isoleucine degradation.


SESI-HRMS data analysis
Data acquisition was completed over a time-span of two years (2020 -2022). The data preprocessing sections are different among experiments because the data preprocessing pipeline was continuously improved.

Data preprocessing -Study I and Study II
Data preprocessing of studies I and II was conducted according to our patented data processing pipeline (European patent No. 20186274.5 and 21185400.5) using MATLAB (version 2021b, Mathworks Inc., USA). To summarize, RAW files were accessed via in-house C# console apps based on Thermo Fisher Scientific's RawFileReader (version 5.0.0.38). To generate the feature list for each experiment, the MATLAB function ksdensity was used for binning the centroid peak list from each scan of all files. XIC were extracted for all features and TIC from each measurement file. For subsequent analysis, only peaks with a signal intensity above 10e3 a.u.
in Thermo's signal intensity scale were further considered when isotopic labeling experiments were performed (study II), or peaks with a signal intensity above 10e4 a.u. for experiments performed in study I. In case satellite artefact peaks were present after applying the defined peak cutoff, those m/z regions were excluded from further analysis (i.e., study I: m/z 50-54, m/z 59-60, m/z 97-98 and m/z 118-119; study II: m/z 50-61, m/z 72, m/z 96-97). Time windows when the samples were pervaded with the gas mixture were determined via a z-scored TIC cutoff. After this, the area under the curves (AUC) during the defined time windows when gas supply through the samples was provided were numerically integrated. As a last step, normalization of AUCs by the time window was performed to get the time-normalized area under the curve (nAUC). Molecular formulae (MF) were assigned based on the "seven golden rules" 1 considering the elements C, H, O, N and S and the adduct [M+H] + . MF-assignment according to the "seven golden rules" 1 is based on an algorithm derived from seven heuristic rules which S4 ensures correct elemental compositions and assigns the most likely and chemically correct MF to a mass spectral feature of interest.
Data preprocessing -Study III MATLAB (version R2022a, MathWorks, USA) was used for data preprocessing. Profile mass spectra and raw centroids (intensity cut-off = 50 2 a.u.) were retrieved by using in-house C# console apps based on Thermo Fisher Scientific's RawFileReader (version 5.0.0.38). Afterwards, reference peaks with formulae fulfilling the "seven golden rules" 1 and common contaminants 2 present in at least 50% of the samples were used to recalibrate centroid and profile mass spectra applying an initial tolerance of 5 ppm. Subsequently, an in-house shape-preserving piecewise cubic interpolation algorithm was deployed to estimate the experimental error across the whole m/z-range. Observed outliers assessed for the reference peaks by moving median algorithm, were excluded when interpolation was performed. This was followed by shifting of the centroids and profile peaks according to the obtained mass error. This process was applied three times for the whole m/z-range of all mass spectra to make sure that the mass error falls below 0.5 ppm. To bin the histograms of the recalibrated centroid peak list, kernel density function was applied. To ensure Gaussian probability density functions of +/-1 ppm at full width at half maximum, an iteration on the bandwidth controlling the smoothness of the probability density curve was performed. Lastly, to create the final data matrix, centroids present in the previously mentioned window were used. Only m/z ≥ 100 were used for subsequent analysis to exclude satellite artefact peaks. MF were assigned based on the "seven golden rules" 1  As a result, from preprocessing, a data matrix of total 216 files (18 files/time points per sample) x 4104 mass spectral features in positive mode was obtained. Time traces for each feature were then computed. Next, the time-normalized AUC of each time trace was calculated for each of the three biological replicates before and after SN addition. To identify differences between the four sample groups after SN addition, a one-way analysis of variance (ANOVA) with AUCs was performed, followed by a posthoc multiple comparison using Tukey's significant difference procedure. For further statistical evaluation only significant features (p ≤ 0.05) between the sample groups of SN1 vs. SN1 + DC and SN2+DC vs. SN2 were further considered. As a second criteria, only features with a log2 fold change (log2FC) ≥1 in mean area under the curve (mAUC) after SN vs. before SN addition were considered as relevant. The data matrix of the remaining features including only the 15 time points after SN addition was then 5 th -root transformed and used to perform t-distributed stochastic neighbor embedding (tSNE) (i.e., Figure   2A and Figure S1A). In addition, the whole data matrix was auto-scaled (z-score) and subjected to hierarchical cluster analysis (Ward method; Euclidean distance) and visualized as heat maps (i.e., Figure 2D, Figure S1B and Figure S1C).

Data postprocessing -Study II
As a result of preprocessing, a data matrix of total 78 files (13 files/time points per sample) x 4384 mass spectral features in positive mode was obtained. Afterwards, features were only considered for further analysis if their signal intensity was greater than zero in at least five out of 13 data time points measured for the samples containing Glc 13 C6. 13 C/ 12 C isotopologue pairs were identified by using a mass difference Δ = n × 1.0034u (± 0.0005u), where n =1, 2, . . . 10 and u = unified atomic mass unit. After this the 13 C/ 12 C ratios were computed for all samples at all time points. Feature ratio time traces containing not a number (NaN) values and higher S6 ratios than expected for natural abundance (natural abundance 13 C/ 12 C ratio = 0.01109) in samples containing Glc 12 C6 were dropped (i.e., ratios > 0.1 as some tolerance due to measurement variability was given). Furthermore, only those features for which the average signal intensity ratios of the samples containing Glc 13 C6 were at least twice as high compared to the average signal intensity ratios of the samples containing Glc 12 C6 were further considered. This ensured that those features showing an increasing ratio trend and therefore an incorporation of 13 C were captured. Features where the suggested number of incorporated 13 C was higher than the proposed MF were not further considered. Only those features which showed consistent 13 C-incorporation among all three biological replicates were finally used for visualizing the 13 C/ 12 C ratio time traces (i.e., Figure 3).

Data postprocessing -Study III
A data matrix of total 20 files (4 biological replicates/group) x 2759 mass spectral features in positive mode and x 938 features in negative mode was obtained. 13 C/ 12 C isotopologue pairs were identified by using a mass difference Δ = n × 1.0034u (± 0.0005u), where n =1, 2, . . . 10 and u = unified atomic mass unit. Next, the signal intensity ratios 13 C/ 12 C were computed for all samples in positive and negative mode. Calculated ratios resulting in NaN values were excluded. ANOVA was then performed for positive and negative ratios of the five groups separately, followed by a post hoc multiple comparison using Tukey's significant difference procedure. Only features fulfilling the following criteria were further considered relevant in negative or positive mode: I) G1, G2 and G3 must be significant (p ≤ 0.05) vs. G4 and G5 or II) G1, G2 and G3 must be significant (p ≤ 0.05) vs. G4 or III) G1, G2 and G3 must be significant (p ≤ 0.05) vs. G5. In addition, only features showing a positive difference in the estimated group means of G4 and G5 vs. G1, G2 and G3 were further considered. Furthermore, if for the significant features, the 13 C/ 12 C ratio in G1, G2 and G3 was higher than expected natural abundance (see data postprocessing study II), the features were not further considered. Features S7 where the suggested number of incorporated 13 C was higher than the proposed MF were also not considered for the final ratio plots (Figure 4, Figure S5 and Table S3).
To perform PCA, original positive and negative matrices were merged. For PCA visualization including all features, data were first log10-transformed and afterwards auto-scaled ( Figure   S6). For PCA visualization including only significant features, a one-way ANOVA test was performed on the log10-transformed data matrix. The data matrix containing the significant features (raw p ≤ 0.05) was then auto-scaled and used to perform PCA ( Figure 5A).
Positive and negative data matrices of G2 and G3 were used for statistical analysis and feature selection via MetaboAnalystR (version 5.0; https://www.metaboanalyst.ca). Data were log10 transformed and pareto-scaled (mean-centered and divided by the square root of the standard deviation of each variable). PCA was then used to first visualize the data ( Figure S7). Statistical analysis was performed using a paired t-test. The top 25 features (Table S4) were then used for cluster analysis and illustrated as a heatmap ( Figure 5B). Furthermore, raw p-values after paired t-test and log2FC were used to perform metabolic pathway enrichment analysis ( Figure 5C,  Putative compound assignment of features presented in Figure 3 and Figure 4 Putative compound assignment was done in MetaboAnalyst R (version 5.0) using mmu_kegg organism library considering only protonated and deprotonated species within 5 ppm.

LC-MS/MS analysis for compound identification of SCH3-BTH
At the end of the volatilome measurements performed in study I with bacterial SN, the content of one cell culture flask containing DCs with SN1 was transferred to a Falcon tube and centrifuged for 5 min at 1500 rounds per min (rpm). After centrifugation, 3 mL of the supernatant was transferred to an Eppendorf tube and stored at -80°C until LC-MS analysis. Standards and samples were analyzed by means of a Vanquish UHPLC system equipped with a Q Exactive