Metaproteomics of saliva identifies human protein markers specific for individuals with periodontitis and dental caries compared to orally healthy controls

Study population and sample collection

The study population, clinical examination and collection of saliva samples have been presented in detail (Belstrøm et al., 2016b). In brief, saliva production was induced by chewing on a tasteless paraffin gum, and chewing-stimulated saliva samples were collected from 10 patients with periodontitis, 10 patients with dental caries and 10 orally healthy individuals following a standardized protocol (Kongstad et al., 2013). Immediately after collection saliva samples were divided into four aliquots and stored at −80 °C for further analysis. One aliquot has previously been analyzed by next-generation sequencing (the Human Oral Microbe Identification using Next Generation Sequencing, HOMINGS) (Belstrøm et al., 2016b). All participants signed an informed consent prior to participation, and the study was approved by the regional ethical committee (H-15000856-53175) and reported to the Danish Data Authorization (2015-54-0970).

Sample preparation

The saliva proteome samples were prepared as described in (Jersie-Christensen, Sultan & Olsen, 2016) with a few modifications. Briefly, 1 ml of saliva was mixed with 1.5 ml lysis buffer (9 M Guanidine hydrochloride, 10 mM Chloroacetamide, 5 mM tris(2-carboxyethyl)phosphine in 100 mM Tris pH 8.5) and heated for 10 min (99 °C) followed by 4 min of sonication.

Protein concentration was measured with Bradford protein assay and ranged from 1–2.5 mg/ml. All samples were digested with the same amount of Lysyl Endoproteinase (Wako, Osaka, Japan) in a ratio of 1:100 w/w calculated from the highest concentration for 2 h. Samples were diluted to a final volume of 10 ml with 25 mM Tris pH8 and digested overnight with Trypsin (modified sequencing grade; Sigma) in a 1:100 w/w ratio.

Digestion was quenched by adding 1 ml of 10% trifluoroacetic acid and centrifuged at 2,000 g for 5 min. The resulting soluble peptides in the supernatant were desalted and concentrated on Waters Sep-Pak reversed-phase C₁₈ cartridges (one per sample) and the tryptic peptide mixtures were eluted with 40% acetonitrile (ACN) followed by 60% ACN. Peptide concentrations were determined by NanoDrop (Thermo, Wilmington, DE, USA) measurement.

Mass spectrometry analysis

A total of 1.5 μg peptide mixture from each sample was analyzed by online nano-scale liquid chromatography tandem mass spectrometry (LC-MS/MS) in turn. Peptides were separated on an in-house packed 50 cm capillary column with 1.9 μm Reprosil-Pur C₁₈ beads using an EASY-nLC 1000 system (Thermo Scientific). The column temperature was maintained at 50 °C using an integrated column oven (PRSO-V1; Sonation GmbH, Biberach, Germany). Buffer A consisted of 0.1% Formic acid, and buffer B of 80% ACN, 0.1% Formic acid. The flow rate of the gradient was 200 nl/min and started at 5% buffer B, going to 25% buffer B in 110 min, followed by a 25 min step going to 40% buffer B and continuing to 80% buffer B in 5 min for a 5 min wash and returning to 5% in 5 min and continuing for re-equilibration for 5 min.

The Q Exactive HF instrument (Thermo Scientific, Bremen, Germany) was run in a data dependent acquisition mode using a top 12 Higher-Collisional Dissociation (HCD)-MS/MS method with the following settings. Spray voltage was set to 2 kV, S-lens RF level at 50, and heated capillary at 275 °C. Full scan resolutions were set to 60,000 at m/z 200 and the scan target was 3 × 10⁶ with a maximum fill time of 20 ms. Full-scan MS mass range was set to 300–1,750 and dynamic exclusion to 20 s. Target value for HCD-MS/MS scans was set at 1 × 10⁵ with a resolution of 30,000 and a maximum fill time of 60 ms. Normalized collision energy was set at 28.

Data analysis

All 30 raw LC-MS/MS data files were processed together using MaxQuant version 1.5.0.38 (Cox & Mann, 2008) with default settings and match between runs. The integrated Andromeda peptide search engine and a reversed database approach applying a 1% FDR at both peptide and protein level was used. The data was searched in two iterations analogous to a previously described metaproteomics database search strategy (Jagtap et al., 2013). First, the search space consisted of the full SwissProt protein database (The UniProt Consortium, 2015) and the Human Oral Microbiome database (Chen et al., 2010) (both downloaded August 2014). The resulting search output was then used for reduction of the search space after filtering on different parameters. As a quality control measure, proteins with less than two unique peptides were removed. Furthermore, we required proteins to be detected in at least five out of 30 samples. Accession numbers from the Majority protein IDs column in the proteinGroups.txt were used to retrieve information about Lowest Common Ancestor (LCA) for each protein group entry. To find the LCA of a protein group, accession numbers with the most peptide-associations were selected, mapped to species and their full taxonomic lineage. The lowest taxonomic rank of the intersection of the latter yielded the LCA. All LCA searches resulting in the parvorder Catarrhini (primates) were set to be human. LCAs at taxonomic rank of species and genera, as well as all of their descendants were used to create a new, reduced search space. The latter was used for the second iteration of MS data identification and quantification and all accession numbers within a protein group were used to perform LCA searches. The above functionality was achieved using the Python programming language. Species names from SwissProt and HOMD were mapped to NCBI taxonomic identifiers using UniProt (http://www.uniprot.org/docs/speclist) and NCBI resources (http://www.ncbi.nlm.nih.gov/Taxonomy/TaxIdentifier/tax_identifier.cgi), respectively. Full taxonomic lineages were retrieved from NCBI Taxonomy database dump files (ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/). Taxonomic comparison at genus- and species level was performed using Mann-Whitney U test with Benjamini-Hochberg correction for multiple testing.

Protein intensities based on summed peptide MS signal intensities were quantile normalized using the limma package version 3.24.15 under R version 3.2.2. Only proteins identified with more than one peptide (“razor + unique”) and present in more than five out of the 30 samples were considered for further analysis. The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium via the PRIDE (Vizcaíno et al., 2016) partner repository with the dataset identifier PXD004319.

For comparative analysis of the human protein profiles, the normalized intensity values were log2 transformed. For principal component analysis (PCA), missing values were replaced with the constant value 19, representing the lowest protein intensity value measured. Analysis of significance (ANOVA) between groups was performed with the software package Perseus (http://www.perseus-framework.org). The resulting differentially expressed proteins were clustered using Euclidian distance after scaling the data by subtracting the mean intensity value. All p-values were corrected for multiple comparisons.

Functional annotation of bacterial proteins

Bacterial proteins from HOMD were searched against Hidden Markov Models (HMMs) of bacterial Nested Orthologous Groups (NOGs) from eggNOG (Huerta-Cepas et al., 2016) using HMMscan version 3.1 (http://hmmer.org) (Eddy, 2009). For each protein query the resulting hits were restricted by two criteria. E-values had to be equal or lower than 1e-4 and a maximum overlap of eight amino acids of HMMs was allowed (selecting hits with the lowest e-value). All corresponding NOG-names were used to retrieve Gene Ontology (GO)-terms as well as KEGG pathways from eggNOG.

KEGG pathway enrichment and characterization

To gain insights into differences between the three sample groups, KEGG pathway enrichment was performed using a modified version of AGOtool (Schölz et al., 2015). Individual samples were grouped to sample categories and the three paired combinations used for the enrichment analysis. All bacterial protein groups with an LCA at rank genus or below were selected. Benjamini-Hochberg correction (FDR) of p-values was applied to correct for multiple testing. The FDR was set to 1%. The following additional filter criterion was applied. The fold change had to be equal or higher than 2 or equal or lower than 0.5.

To get a functional overview of the bacterial proteins, we characterized each individual sample group by counting the number of protein groups associated with each KEGG pathway. For visualization purposes (Fig. S2), we selected the most highly associated terms. Within each group the number of associations was converted to percentages, and the most highly associated terms retained, until a cumulative sum of 90% was reached. This reduced the number of KEGG terms from 135 to 50.

General findings

Biomass analyses based on summed protein intensity measures demonstrated that approximately 95 and 3% of the total protein biomass was of human and bacterial origin, respectively (Fig. 1). Food-related proteins and proteins that could not be assigned to kingdom level comprised the rest of the biomass. We identified a total of 35,664 unique peptides from 4,161 proteins, of which 97% of the identified proteins could be assigned as bacterial or human proteins, with almost equal numbers of the two (Fig. 1; Table 1). The protein biomass and numbers between groups showed no significant differences using ANOVA analysis (Fig. S1).

Human protein profiling

Principal component analysis (Fig. 2A) of the human proteins in saliva showed decent separation of samples from patients with periodontitis and dental caries from orally healthy individuals, based on the most decisive component of the dataset, accounting for 17.9% of the variation. The most enriched KEGG pathway in component 1 and 2 was ‘Complement and coagulation cascades’ (Fig. 2B). Component 2 also separated samples from patients with dental caries and periodontitis patients from orally healthy individuals with the component explaining 12.7% of the variation. Two of the most enriched terms in component 2 in the direction of the orally healthy individuals were KEGG pathway ‘Salivary secretion’ and GOBP ‘protein glycosylation.’

From a total of 2,090 identified proteins of human origin, 60 proteins were significantly differentially expressed when performing multiple sample test (ANOVA, p < 0.05). Hierarchical cluster analysis of the proteins nicely separated the three sample groups, although three periodontitis individuals cluster together with the healthy group (Fig. 3). We identified three main protein clusters. Cluster I contains human proteins that are higher expressed in both disease groups compared to controls, and 10 of 20 proteins in this cluster are associated with the GO term ‘innate immune response’ (protein name in purple). Cluster II consist of nine proteins that distinguish the individuals with caries from the other groups. In cluster III the protein intensities in the caries group are higher than the mean, for the orally healthy group it is around the mean and for the individuals with periodontitis lower.

Bacterial protein profiling

Of the 1,946 proteins of bacterial origin identified, approximately 92% and 34% could be assigned to genus and species level, respectively (Table 1). A total of 29 different bacterial genera and 81 species were identified. The five most predominant bacterial genera were Streptococcus, Prevotella, Veillonella, Rothia and Neisseria collectively representing approx. 70% of the total bacterial mass. The five most predominant bacterial species identified were Rothia mucilanginosa, Veillonella atypica, Prevotella histicola, Prevotella melaninogenica and Streptococcus salivarius. Abundances of the 20 most predominant bacterial genera and species are displayed in Figs. 4A and 4B. No statistically significant differences were observed between groups at genus or species level. However, at genus level there is a trend of higher proportion of Veillonella and lower proportion of Haemophilus were associated with dental caries and higher proportions of Fusobacterium, Leptotrichia and Selenomonas and lower proportions of Streptococcus, Rothia and Haemophilus were associated with periodontitis, when compared to orally healthy individuals. The same trend is seen at species level where higher proportion of Veillonella atypica and lower proportion of Haemophilus parainfluenzae were associated with dental caries, and higher proportions of Fusobacterium periodonticum and Leptotrichia wadei and lower proportions of Haemophilus parainfluenzae were associated with periodontitis, when compared to orally healthy individuals. A full list of all bacterial genera and species identified are presented in Table S1.

KEGG pathway enrichment for bacterial proteins

KEGG pathway enrichment analysis of bacterial proteins resulted in no significant differences with the application of the previously mentioned fold-change and FDR filter criteria. The characterization of functional associations of bacterial proteins is shown in Fig. S2.