Variations in protein expression associated with oral cancer

2.1Subject recruitment and sample collection

Thirty OSCC patients (15 males and 15 females) in the early stage with a mean age of 55 were recruited from the Department of Oral and Maxillofacial-Head and Neck Oncology of the Ninth People’s Hospital (Shanghai, China). The patients did not have detectable periodontal inflammation, oral mucosal diseases, visible carious lesions, or any severe systemic disorders (such as immune compromise, diabetes, or genetic diseases). Furthermore, the patients had not taken antibiotics or received treatment for OSCC at least two weeks prior to the sampling and consented to clinical examination and sampling. The present study was approved by the Ethics Committee of Shanghai Ninth People’s Hospital affiliated with the Shanghai Jiao Tong University School of Medicine. Written informed consent was obtained from all subjects. All experiments were performed according to approved guidelines.

The swabs of oral lesions and anatomically matched normal sites were collected according to a well-defined clinical protocol [13]. Subjects were prevented from drinking and eating for at least 2 h before the sampling. All samples were transported to the laboratory on ice within 2 h of collection and were stored at -80∘C before subsequent processing.

2.2Sample preparation

Use liquid nitrogen to grind the sample into powder in mortar. Quickly transfer all powder samples to 1.5 ml centrifuge tubes precooled in liquid nitrogen, and add 200 to each tube μL SDT protein cleavage buffer (4% SDS, 100 mm Tris HCl and 100 mm dithiothreitol (DTT), pH 8.0); The samples were placed in a boiling water bath for 10 minutes, sonicated in an ice bath for 10 minutes (every 7 seconds, 40 W, 3 seconds), placed in a boiling water bath for 5 minutes, and centrifuged at 14000 g for 30 minutes. Collect the supernatant and pass 0.22 μ M ultrafiltration tube filtration. The filtrate is stored at -80∘C.

After thawing the samples at 4∘C, 10 μL of the samples were supplemented with 2 μL of 5 × loading buffer (250 mM Tris-HCl (pH 6.8), 10% (w/v) SDS, 0.5% (w/v) BPB, 50% (v/v) glycerin, and 5% DTT), mixed, and heated at 98∘C for 3 min; the samples were cooled to room temperature, centrifuged at 3,000 g for 1 min, and loaded to perform electrophoresis (100 V for 10 min and 200 V for 30 min). After electrophoresis, the gel was slightly washed with deionized water to remove electrophoresis buffer. Then, the gels were quickly silver-stained. The images were scanned and recorded. The samples were finally investigated by screening of the images.

2.3Protein extraction and peptide sample preparation for mass spectrometry

100 per sample in total μL for FASP digestion; Add 100 mm DTT to the sample, boil it in water for 5 minutes and cool it to room temperature. Add and mix total 200 μL UA buffer (8 M urea and 150 mm Tris HCl, pH 8.5) and transfer the solution to a 30 kDa ultrafiltration tube; The sample was centrifuged at 14000 g for 30 minutes at room temperature; Discard the filtrate and repeat the above steps 3 times. Add 100 in total μL IAA (50 mm IAA in UA); The samples were shaken at 600 rpm for 1 minute, incubated in the dark at 300 rpm at room temperature for 30 minutes, and centrifuged at 14000 g at room temperature for 30 minutes. Add total 100 μL UA buffer, and centrifuge the sample at 14000 g for 30 minutes at room temperature; These steps were repeated three times. Add 100 in total μL 25 mm ABC buffer and centrifuge the sample at 14000 g for 30 minutes at room temperature; These steps are repeated four times. Finally, discard the filtrate and add 40 μL trypsin buffer (at 40 μ Add 2 to l 100 mm ABC buffer μG trypsin) and place the mixture on a thermostatic mixer (300 rpm, 18 h, 37∘C). Collect the filtrate by centrifuging at 14000 g for 30 minutes at room temperature; Replace the collecting pipe and add 40 μL 25 mm ABC buffer; The sample was centrifuged at 14000 g for 30 minutes at room temperature to obtain the filtrate, and the peptide was quantified using OD 280.

2.4Liquid chromatography-tandem mass spectrometry analysis

Solvent A is an aqueous solution of 0.1% formic acid, and solvent B is a solution of 0.1% formic acid in 100% acetonitrile. Use 95% solvent A to balance the thermal scientific analysis column (75 μm × 25 cm, 5 μm, 100 Å, C18). Use the autosampler to load the sample into the thermoscience easy catch column (100 μm × 2 cm, 5 μm, 100 Å, C18), and then separated by chromatographic column. The relevant liquid phase gradient is as follows: 0–132 minutes, the linear gradient of solvent B ranges from 3% to 24%; 132–200 minutes, the linear gradient of liquid solvent B ranges from 24% to 42%; In 200–235 minutes, the linear gradient of liquid solvent B ranges from 42% to 90%; 90% solvent B is maintained at 235–240 minutes. The enzymolysis products were desalted, separated by capillary high performance liquid chromatography, and analyzed by tandem mass spectrometry using a q-exactive mass spectrometer (thermo Finnigan, San Jose, California). Analysis time: 240 minutes. Detection method: positive ion. Scanning range of parent ion: 350–2000 m/z. Resolution of primary mass spectrometry: 70000@m /z 200. AGC target: 3e6. Level 1 maximum it: 50 ms. Number of scanning ranges: 1. Dynamic exclusion: 30.0S. The mass charge ratio of peptides and peptide fragments was obtained according to the following methods: 20 fragment maps were obtained after each full scan (MS2 scan); MS2 activation type: HCD; Isolation window: 1.2; Resolution of secondary mass spectrometry: 17500@m /z 200; MICROSCAN: 1; Secondary maximum it: 100 ms; AGC target: 1E5; Normalized collision energy: 32 EV; Underfill rate: 2%.

2.5Database creation and searching for identification and quantification of the peptides and proteins

To increase the sensitivity of a large data search, the MetaPro-IQ approach was used to identify and quantify the peptides and proteins using the raw data of mass spectrometry. Briefly, the whole oral microbiome gene catalog database was initially downloaded from the Human Oral Microbiome Database (HOME, version: 9.0, http://www.homd.org). Then, Proteome Discoverer (version 2.1) was used to search against this database to generate a “pseudometaproteome” database for each sample. A reduced database containing all possible proteins derived from the peptide spectrum matches was generated and hyphenated using reversed sequences for each sample for a typical target-decoy database search. FDR < 0.01 filtering was used to confidently confirm the lists of identified peptides and proteins. The final protein database was created by combining and deduplicating all resulting protein lists. Finally, all raw files were imported into MaxQuant (version: 1.6.5.0) for peptide and protein identification and quantification based on this database. Search parameters of Maxquant for the protein database were as follows: trypsin, two maximum missed cleavages; oxidation (methionine) and acetyl (protein N-terminus) as variable modifications; carbamidomethylation (cysteine) as fixed modification; LFQ and iBAQ label-free quantification; 6 ppm for main search; revert decoy mode; PSM FDR 0.01; protein FDR 0.01; 2 minimal peptide length; unique + razor peptides for quantification; 20 ppm FTMS MS/MS tolerance; and 12 FTMS top peaks per 100 Da.

2.6Microbial taxonomic analysis

All peptide sequences identified by MaxQuant were submitted to the ‘Metaproteome analysis’ module of the UniPept web application (version: 4.2) for microbial taxonomic analysis with the “Equate I and L” and “Advanced missing cleavage handling” options. Moreover, all identified top rank proteins were assigned to the corresponding taxa by an in-house Perl script based on the taxonomic level file obtained from the HOMD. Unmatched top proteins were manually searched using the Taxonomy database of the NCBI by keywords (genus name) in protein description. Statistics and plots of taxonomic analysis were assessed by R (version: 3.5.1, https://www.r-project.org/). The Krona charts of identified bacteria were created by the Krona Excel template based on protein counts and LFQ intensity accumulation (Krona Excel template: https://github.com/marbl/Krona/wiki/ExcelTemplate).

At all taxa levels, the Mann-Whitney U test and paired sample t test were used to identify differentially abundant taxa in the two groups based on log10-transformed cumulative LFQ intensity calculated using R.

2.7Microbial function and pathway analysis

A standalone BLASTP(version:2.2.31+) was used to match all identified microbial proteins to the UniProtKB/Swiss-Prot database, which included 560,118 sequences, and the COG database (version: 2014), which included 1,781,653 sequences for obtaining function annotations. Furthermore, gene ontology (GO) functional annotations were predicted by InterProScan (version: 5.35–74.0) with default parameters. Moreover, KEGG Orthology Based Annotation System (KOBAS version:3.0), parameters: E-value<=1e-05, rank<=5) was performed for the KEGG Orthology (KO) annotation prediction and binomial tests were performed between the B and Z groups at the KEGG pathway level. Then, the KronaTools was used to plot the Krona chart for visualizing assignment of the proteins at the KEGG pathway level. Besides, statistical analysis and graphics drawing were completed with R (version:3.6.3) .

2.8Statistical and multivariate data analysis

We used the equation log10(LFQ intensity + 1) to transform the LFQ intensity of quantified protein groups in R for subsequent statistical analysis. After removing the proteins that lacked LFQ signals in all samples, univariate significant differences between the B and Z groups were assessed by the paired sample T test using R (version: 3.5.1). Moreover, after imputing missing values as one thousandth of the minimum value, partial least squares discriminant analysis (PLS-DA) was performed to assess differentially expressed proteins using the MetaboAnalystR package. Only those protein groups with valid LFQ intensity values in ⩾ 50% of the samples were used for multivariate statistical analysis. Principal component analysis (PCA) and PLS-DA were performed by the mixOmics package. Hierarchical clustering analysis and heatmap plotting were performed by the R pheatmap package.

2.9Verification with parallel reaction monitoring technology

The appropriated target peptides were selected by analyzing the results of the original label-free-based quantitative. Then, these peptides were performed targeted shotgun MS to select peptides for PRM quantification analysis. Finally, the Xcalibur software was used to implement Parallel Reaction Monitoring (PRM) analysis of target peptides.

Statistical analysis was performed by R and the ‘MetaboAnalystR’ package.

2.10Data availability

All MS raw files were submitted to the ProteomeXchange Consortium (http://www.proteomexchange. com) via the PRIDE partner repository with the dataset identifier.

3.1Statistics of the number of identified peptide sequences and protein groups

We quantified a total of 38,536 peptide sequences and 11,469 protein groups that covered at least one peptide in 60 samples of 30 OSCC patients. A total of 26,389 (68.48%) of these peptide sequences were confirmed as unique peptides. On average, 1,428 ± 304 peptide sequences and 1,382 ± 194 protein groups were identified per sample. In addition, a mean of 1,790 ± 503 peptide sequences and 1,632 ± 289 protein groups were identified per sample of Group B (swabs of lesion surfaces). As expected, 1,066 ± 336 mean peptide sequences and 1,131 ± 256 protein groups for Group Z (contralateral normal mucosa) were slightly lower than the counts detected in Group B (Fig. 1). Using the ‘Stats’ package of R, p= 0.009108 and p value = 0.03148 was computed based on the paired sample T test of the B versus Z groups. The p values indicated that the peptide and protein counts of Group B were significantly greater than those of Group Z. A total of 10,773 nonredundant protein groups were identified in Group B samples, and 7,855 nonredundant protein groups were identified in Group Z samples and had 7,159 shared proteins between the two groups.

Figure 1.

Peptide and protein group identification by oral metaproteomics. (A) Box plot representations of the number of peptides identified in the present study. (B) Box plot representations of the number of identified proteins. Different colors of the symbols are used to represent individual samples, and the median (central lines), mean (white small dot), 25% and 75% quartile ranges (box width), and upper and lower limits (asterisk) are shown. All: All samples from 30 OSCC patients; B: Swabs of lesion surfaces; Z: Contralateral normal mucosa.

3.2Taxonomic analysis of the microbiome by metaproteomics

Taxonomic characterization is important in microbiome studies. In the present study, we attempted to use two strategies to explore the taxonomic composition of the microbiome. According to the first strategy, all 38,536 identified peptide sequences were submitted to UniPept (version: 4.2) to perform the lowest common ancestor (LCA) analysis. A total of 37,475 (98.79%) peptides were matched by the UniPept. A sunburst diagram was generated to provide a comprehensive view of the four taxonomic levels (Fig. 2).

Figure 2.

Sunburst diagram for bacteria taxonomic composition.

In addition to taxonomic analysis, matched peptides were related to enzymes, and gene ontology (GO) functional categories included GO biological processes (GOBP), GO cellular components (GOCC), and GO molecular functions (GOMF) in UniPept. A total of 1,326 enzymes, 966 GOBP terms, 255 GOCC terms, and 1,288 GOMF terms were matched. The top 20 and top 10 peptide count-related enzymes and GO categories are shown as bar charts (Fig. 3).

Figure 3A.

Top 20 of peptide counts related to enzymes.

Figure 3B.

Top 10 peptide counts related to GOBP terms.

Figure 3C.

Top 10 peptide counts related to GOCC terms.

Figure 3D.

Top 10 peptide counts related to GOMF terms.

The second strategy involved assignment of all identified top proteins assigned to the taxa using an in-house Perl script to search for top protein descriptions using the taxonomic level file downloaded from the HOMD. A total of 10 phyla, 19 classes, 33 orders, 60 families, 113 genera, and 261 species were identified in the present study based on at least one top protein (Fig. 4).

Figure 4.

Protein count of taxonomic levels.

Using cumulative LFQ intensity transformed by log10, the Mann-Whitney U test and paired sample T test were performed at the phylum, class, order, family, genus, and species levels. Based on the values of the Mann-Whitney U test of p< 0.05 or of paired sample T test of p< 0.05, 68 significantly different taxa were identified, including 2 phyla, 3 classes, 6 orders, 11 families, 12 genera, and 34 species. These taxa included Eubacteriaceae [XV] at the family level and Selenomonas flueggei at the species level, which were identified only by the Mann-Whitney U test, and 38 taxa were identified only by the paired sample T test method; 28 taxa were identified by the two methods (Table 1).

3.3COG functional analysis of identified microbial proteins

A total of 10,696 (93.26%) identified microbial proteins were annotated by using BLASTP to search against the COG database, and 1,388 COGs were obtained based on only the best hit for each protein. These proteins were assigned to 23 COG categories in five top level classifications, including 1,386 proteins in “Energy production and conversion”, 146 proteins in “Cell cycle control, cell division, chromosome partitioning”, 1,098 proteins in “Amino acid transport and metabolism”, 287 proteins in “Nucleotide transport and metabolism”, 1,335 proteins in “Carbohydrate transport and metabolism”, 259 proteins in “Coenzyme transport and metabolism”, 408 proteins in “Lipid transport and metabolism”, 2,364 proteins in “Translation, ribosomal structure and biogenesis”, 303 proteins in “Transcription”, 184 proteins in “Replication, recombination and repair”, 637 proteins in “Posttranslational modification, protein turnover, chaperones”, 113 proteins in “Cell motility”, 698 proteins in “Posttranslational modification, protein turnover, chaperones”, 484 proteins in “Inorganic ion transport and metabolism”, 74 proteins in “Secondary metabolites biosynthesis, transport and catabolism”, 294 proteins in “General function prediction only”, 210 proteins in “Function unknown”, 105 proteins in “Signal transduction mechanisms”, 177 proteins in “Intracellular trafficking, secretion, and vesicular transport”, 122 proteins in “Defense mechanisms”, 2 proteins in “Extracellular structures”, 6 proteins in “Mobilome: prophages, transposons”, and 4 proteins in “Cytoskeleton” (Fig. 5).

Figure 5.

Bar chart for COG functional classification.

At the COG category level, we used the same methods to identify significant differentially abundant categories. Unexpectedly, no significant differentially abundant categories were detected. All p values of the Mann-Whitney U test were greater than 0.05, and the maximum VIP value was less than 2 (R: general function prediction only; VIP = 1.9795).

Figure 6.

Pie charts of protein count at the GO level 2 of biological process, cellular component, and molecular function ontology.

3.4Gene ontology functional analysis

Overall, 9194 (80.16%) top identified microbial proteins in general, according to interproscan and internal Perl scripts, these top identified microbial proteins correspond to at least one go term. The number of proteins is counted at the GO level 2 of biological processes, cellular components and molecular functional Ontology (Fig. 6). Then, Fisher’s exact test was performed to compare the difference in the number of proteomes between group B and group Z. Based on p< 0.05, GO terms in group B were significantly richer than those in group Z, including 35 entries in biological process ontology, 2 entries in cell component ontology, and 30 entries in molecular function ontology. GO terms in group Z were significantly richer than those in group B, including 59 entries in biological process, 17 entries in cell component, and 30 entries in molecular function ontology (Fig. 7).

Figure 7.

Bubble charts of top 10 GO terms. Left: Top 10 GO terms of Group B significantly more abundant than those of Group Z. Right: Top 10 GO terms of Group Z significantly more abundant than those of Group B. Y-axis corresponds to significant GO terms, and X-axis corresponds to negative log10 transformed p value according to the results of the Fisher’s exact test. The dot size corresponds to protein counts.

Use internal Perl and R scripts to enumerate the proteins corresponding to the go items of all samples. Fisher’s exact test was performed using R to determine important GO items based on the number of proteins in paired samples from the same patient. At p< 0.05, in at least 12 paired samples, 80 GO terms in group Z (33 terms in biological process ontology, 5 terms in cell component ontology, 42 terms in molecular function ontology) were significantly richer than those in group B. compared with group Z, only 10 GO terms in group B were richer (8 in biological process ontology, 2 in molecular function ontology). Then, we use the negative logarithm 10 to convert the p value and use the heat map generated by R to visualize the data (Fig. 8).

Figure 8.

Heatmaps of transformed p values of significant GO terms in at least 12 sample pairs. A) Significant GO terms of the Z group were more abundant than those of the B group; B) Significant GO terms of the B group were more abundant than those of the Z group.

3.5Pathway analysis based on KEGG

A total of 9,039 top identified microbial proteins were annotated in the KEGG database using KOBAS 3.0. These proteins were related to 1,385 nonredundant KO identifiers and assigned to 196 KEGG pathways belonging to 41 KEGG categories of the 6 top classifications (Fig. 9).

Figure 9.

Krona plot of proteins corresponding to the KEGG pathways. Krona plot shows the proportion of proteins assigned to the KEGG pathways from the top levels to specific pathways.

Significantly different pathways with p< 0.05, according to the binomial test via KOBAS 3.0 based on the protein counts in the pathways, were identified as 8 pathways of Group B more abundant than those in Group Z. In these pathways, the FDR-corrected p values of the bacterial secretion system (FDR = 0.017) and protein export (FDR = 0.024) were less than 0.05. Moreover, we identified 18 significantly different pathways of Group Z more abundant than those of Group B using the same method. Nine pathways had an FDR < 0.05, including glycolysis/gluconeogenesis (FDR = 0.0016), ribosomes (FDR = 0.0058), carbon metabolism (FDR = 0.0058), RNA degradation (FDR = 0.0058), biosynthesis of amino acids (FDR = 0.0108), HIF-1 signaling pathway (0.0108), methane metabolism (FDR = 0.0138), tuberculosis (FDR = 0.0189), and carbon fixation in photosynthetic organisms (FDR = 0.0232). Moreover, a paired sample test used the same methods to obtain significantly different pathways. A total of 72 pathways with at least one p value of paired sample test less than 0.05 were identified in the case of Group B vs. Group Z, and a total of 88 pathways were identified in the case of Group Z vs. Group B. Some pathways were significantly more abundant in Group B than those in Group Z based on the binomial test selection of at least 9 paired sample tests with a p value < 0.05, including carbon fixation in photosynthetic organisms (9), HIF-1 signaling pathway (9), ribosome (9), and carbon metabolism (10); carbon fixation in photosynthetic organisms (9), oxidative phosphorylation (9), Alzheimer’s disease (10), methane metabolism (10), biosynthesis of amino acids(11), glycolysis/gluconeogenesis (12), HIF-1 signaling pathway (12), photosynthesis (12), longevity regulating pathway – worm (14), tuberculosis (14), and RNA degradation (17) were more abundant in Group Z than the corresponding pathways in Group B. (Table 2).

Figure 10A.

Venn diagram for candidate DEPs identified by PLS-DA (VIP), Mann-Whitney U test (MWT), and paired sample T test (PSTT).

Figure 10B.

Heatmap of DEPs. Row clustering includes Manhattan distance and complete linkage.

3.6Selection of differentially expressed proteins

Log10-transformed LFQ intensities were used to identify a total of 39 protein groups as candidate significantly differentially expressed proteins (DEPs) in the B versus Z groups according to the Mann-Whitney U test at a cutoff p value < 0.05. Moreover, a total of 93 candidate DEPs were identified using a paired sample T test with a cutoff p value < 0.05. Furthermore, missing values were imputed as one thousandth of the minimum values, and PLS-DA was performed. A total of 399 candidate DEPs were assessed based on VIP > 2. DEPs obtained by the three methods were used to plot a Venn diagram (Fig. 10A), and 38 common DEPs were finally identified and illustrated on a heatmap generated by the pheatmap package (Fig. 10B). Receiver operating characteristic (ROC) analysis was performed on 38 common candidate DEPs using the MetaboAnalystR package.

Figure 11.

PCA score plot of proteins quantified in ⩾ 50% of the samples. B) PLS-DA score plot of proteins quantified in ⩾ 50% of the samples.

3.7Multidimensional variable analysis

We selected proteins quantified in ⩾ 50% of the samples for principal component analysis (PCA) and partial least squares discriminant analysis (PLS-DA). A total of 10 proteins were selected. PCA and PLS-DA were performed using these proteins by the mixOmics package in R. (Fig. 11).

3.8Mass spectrometry analysis of PRM validation candidate peptides and proteins

The results of PRM indicated that the 22Z and 22B samples were unsuitable for analysis. Therefore, these two samples were removed from the analysis. Finally, the results of 56 samples from 28 patients were used to perform differentially expressed peptides and protein analysis. We used R (version 3.5.1) software to perform paired sample T test for the selection of differentially expressed peptides and proteins, respectively.

The results of the analysis of proteins identified 10 and 11 candidate differentially expressed proteins according to the paired sample T test with p< 0.05, respectively. The combination of candidate differentially expressed proteins included 2 proteins (hinf_c_1_1270 and fnuc2539_c_1_361) that were determined to be candidate differentially expressed proteins according to metaproteomics analysis [17] (Table 3), and the intersection included 5 proteins.