Metabolomics of pulmonary exacerbations reveals the personalized nature of cystic fibrosis disease

Sample collection

Sputum samples were collected at the adult cystic fibrosis center at the University of California at San Diego Medical Center as described in Lim et al. (2014). Procedures were approved under the UCSD Human Research Protections Program protocol #081500. Treating physicians determined which patients would be selected for the study, whether or not they were presenting with an exacerbation (principally based on the Fuchs criteria (Fuchs et al., 1994)), and the subsequent treatment regime. Samples were classified as ‘exacerbation,’ according to the clinicians diagnosis of the need for intravenous antibiotic administration, ‘treatment,’ as any sample collected during the 14 day treatment course, ‘post-treatment,’ as samples collected immediately after treatment course, and ‘stable,’ as samples collected during routine clinical visits that were not known at the time to be during an exacerbation. Sample collection involved an initial saline mouth rinse followed by expectoration of sputum into a sputum cup after inhalation of 7% hypertonic saline for 30 min and then syringe homogenized according to Lim et al. (2013). Two different sample sets were collected for assessment of metabolome similarity in different clinical states or between patients, one in 2012 and another in 2014 (Table S1). For a separate and more in depth longitudinal analysis, another sample set was collected from a single patient through 1,492 days (n = 37 sputa, Table S2). This sample set captured four separate exacerbation events, including samples that were collected daily for 14 days through treatment of the second exacerbation. Bacterial 16S rRNA gene profiles were previously published on sputum samples from this daily collection of the second CFPE event (Quinn et al., 2015a) and this data is also utilized in this study (see methods below and Table S2).

Extraction and LC-MS/MS

The samples were thawed and then extracted using a sequential method of ethyl acetate solvation, followed by methanol solvation. A volume of 200 µl of sputum was first mixed with 200 µl of ethyl acetate, extracted for 1 h at room temperature and then briefly centrifuged. The supernatant was decanted and then evaporated in a centrifugal evaporator. The same volume of methanol was then added to the remaining sputum sample and incubated at room temperature for 1 h, and then briefly centrifuged. This supernatant was added to the dried pellet of the ethyl acetate extract and then evaporated in a centrifugal evaporator. All pellets were solubilized for 1 h in 100 µl of methanol prior to LC-MS/MS analysis. Mass spectrometry was performed using a Bruker Daltonics^® Maxis qTOF mass spectrometer equipped with a standard electrospray ionization source. A water/acetonitrile solvent separation gradient containing 0.1% formic acid was used from 98:2 to 2:98 water:acetonitrile for a total run time of 840 s. The mass spectrometer was operated in data dependent positive ion mode, automatically switching between full scan MS and MS/MS acquisitions.

Metabolome generation and statistical analysis

For identification of the number of unique spectra in each clinical disease state MS/MS spectral alignments using molecular networking on GNPS (gnps.ucsd.edu) was performed on data from the 2012 and 2014 longitudinal collections. The number of unique spectra identified from molecular networking for each disease state was calculated and visualized using a Venn diagram.

Statistical analysis including multivariate comparisons and quantification of metabolite relative abundance was completed using MS¹ traces. Molecular features detected in the mass spectrometer (MS¹ level) were identified using the Bruker Daltonics^® DataAnalysis software and imported into the ProfileAnalyisis software version 2.1 (build 282, 64-bit) for metabolome generation on the entire data set. The 2012 and 2014 sputum datasets were processed and analyzed separately due to batch effects, and global similarity between patient and clinical state was only compared within each dataset. The metabolomes for dataset 1 and 2 were analyzed with a Bray–Curtis distance matrix separately using the vegan package in R (Oksanen et al., 2015). Each distance matrix was then projected with nMDS to visualize clinical state and patient metabolome similarity. The distance matrices were tested for similarity between the two classifiers using ANOSIM and PerMANOVA (adonis) from the vegan package in the R statistical software. The ANOSIM R statistic is a permutation test of the null hypothesis that within group variation is not greater than between group variation; an R-value above 0.4 is considered sufficient to reject the null hypothesis. The PerMANOVA test is a non-parametric multivariate analysis of variance using a pseudo-F-test. Here the larger the value of F the greater likelihood the null hypothesis of no differences between group variations is false and the p-value comes from permutations. We used the ANOSIM R and PerMANOVA F to test whether the metabolomic data classified more strongly by patient or clinical state.

To identify global biomarkers of CFPE, it was initially necessary to unify MS¹ data from all cohorts and batches into a contiguous matrix containing spectral features from all samples. In order to do so a novel algorithm was used, which employs robust methodologies to overcome difficulties in identifying the same metabolite across multiple data sets (details in Supplemental Information 1). With this unified MS¹ spectral data matrix, a supervised random forest using clinical state as the classifier produced a variable importance plot that identified differentially abundant features across clinical states. Ten features were identified as especially able to differentiate between the clinical states. These features were tested for significant differences using an ANOVA and Tukey’s test of significance (p < 0.05) and corrected for multiple comparisons. Biomarkers in the CF1 longitudinal dataset were identified first using a random forests classification based on disease state. The variable importance plot (VIP) of these random forests was then used to identify the 30 most differential features across the clinical states. These 30 features were then subjected to a one-way ANOVA and a Tukey’s test of significance with a Bonferroni correction (30 features compared for CF1 dataset, p < 0.00167).

16S rRNA gene sequencing and metabolome correlations

Bacterial relative abundance data at the level of genus that was previously generated on the same samples collected during the second exacerbation of the CF1 longitudinal collection (published in reference Quinn et al., 2015a) was used to compare correlations between the metabolome and microbiome variables. A correlation matrix was calculated such that the area under the curve for the abundance of each molecular feature was regressed against the normalized abundance of either Pseudomonas, Rothia, Streptococcus, Stenotrophomonas or total anaerobes (Prevotella, Veillonella, Gamella and Oribacterium). The Pearson correlation coefficient (r) of this regression was then used to build a network in the Cytoscape^® software.

Metabolome annotation

LC-MS/MS files generated on the mass spectrometer were converted to the .mzXML format for dereplication and molecular networking. Spectra generated from the MS/MS acquisition were searched against the GNPS database (gnps.ucsd.edu) using the molecular networking algorithm (Watrous et al., 2012).

Metabolome relationships between patient and clinical state

Two sputum sample sets were collected for this study, one in 2011–2012 (n = 25 sputa, 6 patients) and another in 2014 (19 sputa, 5 patients, Table S1). Both data sets contain at least three samples per patient collected during one of 4 clinical states (Exacerbation ‘Ex,’ Treatment ‘Tr,’ Post Treatment ‘Pt,’ and Stable ‘St’). Technical variability of this method within a single run was assessed using multiple extractions of the same sample and on different days, revealing a highly reproducible method within a LC-MS/MS batch (Fig. S1, Supplemental Information 1 and 3). The data was analyzed with the molecular networking algorithm that identifies unique MS/MS spectra (a proxy for molecules) available through the Global Natural Product Social Molecular Networking database (GNPS, gnps.ucsd.edu, (Watrous et al., 2012)). A total of 4,639 unique MS/MS spectra were detected in the sample set. Exacerbation samples had 556 unique metabolites, Tr samples contained 132, Pt samples contained 781, St samples contained 100, and 1,222 metabolites were common across all clinical states (Fig. 1A).

Batch effects were detected between the sample sets by clustering of the metabolome data (Fig. S1). Due to these LC-MS/MS batch effects, the similarity of each metabolome could not be compared across the two datasets, only within each dataset. Therefore, a Bray–Curtis distance matrix was generated on the metabolite abundance matrix for each sample set separately and visualized with non-metric multidimensional scaling (nMDS) (Fig. 1B). The metabolomes were tested for similarity within the patient or clinical state classifiers using the analysis of similarity (ANOSIM) and permutation multivariate analysis of variance (PerMANOVA) (Fig. 1B and Table 1). The metabolomes clustered more strongly by patient than clinical state for both datasets (Fig. 1B). The ANOSIM statistic R (a measure of the difference of mean ranks between and within groups) was higher by patient classification than clinical state (for both datasets patient R > 0.4, Table 1). The PerMANOVA F-test verified this trend (for both datasets patient F ≥ 2.0, Table 1). This untargeted metabolomics analysis demonstrated that although there were a large number of unique molecules belonging to each clinical state, particularly exacerbation, the overall data did not reveal similarity in these clinical states across patients.

Global biomarkers of CFPE

Common features across the two sample sets were identified using an advanced retention time alignment and feature detection method (Fig. S1, Supplemental Information 1 and 3). Although the majority of metabolites unique to exacerbation were patient specific, some were significantly elevated across patients. These were identified using a supervised random forests (Liaw & Wiener, 2002) approach with the samples classified by clinical state on the merged abundance matrix between the two sample sets. Data from both sets were merged using a novel algorithm expressly designed to overcome inter-batch m∕z and retention time variations in order to maximize the number of metabolites common to both sets (see Supplemental Information 1). The variable importance plot revealed eight metabolites had a mean decrease in accuracy of classification above 4, indicating they were highly enriched in a particular disease state, the remaining metabolites contributed less to the classification (Fig. S2). Of these eight metabolites, MS/MS spectral matching and GNPS library searching allowed putative annotations for two. Platelet activating factor (m∕z524.36 (M + H⁺)) was significantly more abundant in Ex than Tr and Pt states, but not during St (Tukey’s test of ANOVA, Ex–Tr p = 0.021, Ex–Pt p = 0.022, Fig. 2). The diacylglycerophosphocholine lipid PC(18:0/3:1) was significantly more abundant in Ex than all other clinical states (Ex–Tr p = 0.0062, Ex–Pt p = 0.003, Ex–St p = 0.038) (Fig. 2). Lyso-PAF, a related metabolite, was not significantly different between clinical states.

Personalized biomarkers of CFPE

In light of the individuality observed between CF sputum metabolomes, the remainder of this study focused on a single patient with longitudinal sampling over a four-year time period, which captured four separate exacerbation events (CF1 dataset, Table S2, n = 37 sputa). Known spectra detected in the CF1 longitudinal dataset were identified through the GNPS molecular networking workflow and library searching and their abundances tracked through time (Figs. 3A and 3B). The abundance of PAF, Lyso-PAF and a monoacylglycerophosphocholine lipid Lyso-PC (1-O-hexadecyl-2-C-methyl-3-phosphatidylcholine) were monitored through the data set and revealed fluctuations in the abundance of these metabolites through time (Fig. 3A). PAF and Lyso-PAF were elevated during Ex and Tr and the Lyso-PC was elevated during a 6-month period of stability (Fig. 3A). A Tukey’s test of a one-way ANOVA for each molecule across the clinical states Ex, Tr and St (no Pt samples were available from this patient) showed that PAF and Lyso-PAF were significantly more abundant during Ex than during St (PAF p = 0.035, Lyso-PAF p = 0.039), but not Tr (Fig. 3C). PAF and Lyso-PAF abundance was highly correlated through the sampling period (Pearson’s r = 0.96, p < 0.0001), but not Lyso-PC (Pearson’s r = 0.25, p = 0.13). Sphingolipids were also identified in this patient through GNPS library searching including ceramide, sphingomyelin and sphingosine. Ceramide was especially elevated during treatment of an exacerbation (Fig. 3B). The Tukey’s test revealed that this molecule was significantly more abundant during Tr compared to St (p = 0.0008, Fig. 3D). Ceramide and sphingosine were significantly correlated in their abundances (Pearson’s r = 0.53, p = 0.0007), but sphingomyelin and ceramide were not (Pearson’s r = 0.27, p = 0.10). These known molecules may represent specific biomarkers of CF1 associated with inflammation during the onset of an exacerbation and through its treatment.

Other molecules that were not annotated through GNPS were also statistically significant biomarkers of CFPE for this patient. Thirty of these metabolites were identified using a random forests based classification of the entire CF1 longitudinal data set based on clinical state (Fig. S3). Molecular features with parent masses m∕z441.317, m∕z304.280, m∕z508.331, and m∕z551.354 were all significantly elevated during Ex and Tr, compared to St, after a Bonferroni correction for a Tukey’s test of a one-way ANOVA (p < 0.00167, 30 comparisons, Figs. S4 and S5). Metabolites m∕z508.331 and m∕z551.354 both contained a phosphocholine fragment in their MS/MS spectra (Fig. S4), indicating that they are putatively a monoacylglycerophosphocholine and a glycerophospholipid of unknown structure, respectively. Although these metabolites remain unidentified these molecules are especially strong biomarkers of exacerbation in CF1.

Xenobiotic dynamics in CF1

The antibiotics azithromycin, trimethoprim, sulfamethoxazole, and linezolid, and the bronchodilator albuterol were detected in the longitudinal sampling data of CF1. There was varying abundance of these antibiotics that fluctuated through time (Fig. 4). Azithromycin, trimethoprim, sulfamethoxazole, and albuterol were administered consistently as ‘suppressive’ therapy to this patient throughout the sampling period. During the second exacerbation event, trimethoprim, sulfamethoxazole and linezolid, were administered intravenously. The abundance of azithromycin was high in sputum through most of the sampling period, but fluctuated during suppressive administration, occasionally out of the detectable range. As expected, the abundance of trimethroprim and linezolid increased during intravenous therapy, but sulfamethoxazole remained at a very low abundance. Albuterol also increased in abundance during the second exacerbation (this metabolite was not part of the intravenous treatment regime, but was administered), while being relatively low during chronic oral therapy.

Microbial and metabolome correlations

The second exacerbation event during the CF1 longitudinal sampling also had bacterial 16S rDNA amplicon sequencing data previously published in a separate manuscript on the same sputum samples (Quinn et al., 2015a). Therefore, abundances of Pseudomonas, Stenotrophomonas, Streptococcus, Rothia and anaerobes (Veillonella, Prevotella, Gemella and Oribacterium) from the published data were regressed against the metabolite abundances in the same CFPE1 samples (n = 16). Statistically significant Pearson correlations (p < 0.05) were visualized by networking analysis using Pearson’s r as the edge values, such that molecules that were positively correlated with a particular bacterial genus would be connected in the network (Fig. 5). The network topology revealed that there were 168 molecular features significantly correlated with Pseudomonas abundance, 66 with Stenotrophomonas, 46 with Streptococcus, 23 with Rothia, and 22 with anaerobes. The network topology showed two groups of molecules, those that correlated with either Stenotrophomonas or Pseudomonas and those that correlated with Streptococcus, Rothia, and anaerobes (Fig. 5). A number of molecules were correlated to more than one bacterium, but there were no correlations between these two groups of bacteria. This indicated that there were two microbial and molecular communities in this patient during this exacerbation, one associated with Stenotrophomonas and Pseudomonas and one with Gram-positives and anaerobes. These two communities have been previously identified as the Climax and Attack communities, respectively, from a parallel study of these samples (Quinn et al., 2015a).

Automatic metabolite annotation through GNPS led to putative identification of 14 of the 303 molecules in the correlation network. Ceramide (18:2/16:0) and sphingomyelin (16:1/16:0) were correlated with Streptococcus, Rothia and anaerobes. N,N-dimethyl arachidonly amide and arachidonyl amide were correlated to Rothia and anaerobes, but not Streptococcus. Two diacylglycerophosphocholines were correlated with presence of Streptococcus. A number of antibiotics were correlated with Stenotrophomonas and Pseudomonas including trimethoprim, linezolid and acetylsulfamethoxazole, all three were administered intravenously as treatment for this exacerbation. Cholesterol, tryptophan, another diacylglycerophosphocholine, and ceramide (18:1/24:0) were correlated specifically with Pseudomonas relative abundance.