Proteomics and metabolomics studies in pigmented villonodular synovitis uncover the regulation of monocyte differentiation by the ADGRE5-NF-κB pathway | BMC Medicine

Patients and clinical samples

This study was approved by the Biomedical Research Ethics Committee of Sichuan University (No. 125 2020-(921)). The baseline characteristics of the patients are provided in Table 1. Arthroscopic images showed red-brown synovial tissue in PVNS patients, which was markedly hyperplastic compared to the control group due to hemosiderin deposition (Fig. 1A). MRI revealed that the synovial tissue in PVNS patients displayed high signal intensity in the fat-saturated proton density sequence (Fig. 1B). The knee joints of PVNS patients exhibited severe joint swelling and cartilage erosion (Fig. 1B).

Biological samples: We enrolled ten patients with knee PVNS who were admitted to our hospital between December 1, 2019, and May 30, 2020, and met the inclusion criteria. We retrospectively analyzed these patients’ clinical characteristics, imaging changes, and arthroscopic findings. The PVNS group specimens were obtained from patients diagnosed with knee PVNS through postoperative pathological examination, while the control group specimens were obtained from patients with meniscus or ligament injuries. Patients with infections, autoimmune diseases, and metabolic disorders were excluded. Synovial fluid from both groups was collected during arthroscopy, centrifuged, and stored in liquid nitrogen, then transported and frozen at − 80 °C for future use. General clinical and demographic data for the PVNS group (P) and the relative control group (C) are shown in Table 1. Patients with PVNS undergo arthroscopic synovectomy to remove the affected synovial tissue. For meniscus injuries, the treatment is based on the severity of the damage, with options including meniscus repair or meniscectomy. Ligament injuries are treated according to the severity, with options for ligament repair or reconstruction surgery. Among the 10 PVNS patients (mean age = 36.10 ± 11.26 years), 2 (20%) were male. The mean body mass index (BMI) (23.86 ± 2.885 vs. 22.60 ± 2.382 kg/m², p > 0.05), C-reactive protein (CRP) (mg/L, 5.751 ± 3.947 vs. 1.915 ± 1.053; p = 0.0137), triglycerides (TG) (mmol/L, 2.216 ± 1.303 vs. 1.151 ± 0.7565; p = 0.0383), and erythrocyte sedimentation rate (ESR) (28.10 ± 18.42 vs. 8.400 ± 5.481; p = 0.0082) were higher in the PVNS group than in the control group. Serum biochemical parameters, including white blood cell count (WBC), high/low-density lipoprotein (HDL/LDL), total cholesterol (TC), absolute monocyte count (MO#), and absolute lymphocyte count (LY#), were measured using an automated serum biochemical analyser.

Proteomics

Sample preparation and protein extraction: synovial fluid samples (40 μL each) were mixed with 250 μL phosphate-buffered saline (PBS) containing a protease inhibitor (Roche, 4,693,132,001) to prevent protein degradation. Proteins were denatured by adding 250 μL of ST buffer (2% SDS, 100 mmol/L Tris–HCl, pH 7.6), followed by centrifugation at 8000 × g for 1 min. The supernatant was collected, boiled for 5 min, and sonicated to ensure complete protein solubilization. After a second centrifugation (8000 × g, 15 min), the supernatant was collected, and protein concentration was determined using a bicinchoninic acid (BCA) assay.

Protein digestion and peptide preparation: proteins were reduced with 10 mmol/L dithiothreitol (DTT) at 56 °C for 1 h and alkylated with 20 mmol/L iodoacetamide (IAA) in the dark for 30 min at room temperature. The solution was buffer-exchanged into 50 mmol/L ammonium bicarbonate using a FASP column (PALL, OD010C34) and digested overnight with trypsin at 37 °C (enzyme-to-protein ratio, 1:50). The resulting peptides were quantified using a fluorometric peptide assay (Thermo Fisher Scientific, 23,275).

High-pH Reverse Phase Fractionation: Equal amounts of peptides from all enzymatically digested samples were pooled and fractionated using an Agilent 1100 HPLC system with a mobile phase at pH 10. Separation was performed using an Agilent Zorbax Extend-C18 column (2.1 × 150 mm, 5 μm) with UV detection at 210 nm and 280 nm. The mobile phases consisted of the following: mobile phase A: ACN-H₂O (2:98, v/v); mobile phase B: ACN-H₂O (90:10, v/v); both mobile phases were adjusted to pH 10 with ammonium hydroxide.

The gradient elution program was as follows: 0–10 min: 2% B (isocratic); 10–10.01 min: 2–5% B; 10.01–37 min: 5–20% B (linear gradient); 37–48 min: 20–40% B (linear gradient); 48–48.01 min: 40–90% B; 48.01–58 min: 90% B (isocratic); 58–58.01 min: 90–2% B; 58.01–63 min: 2% B (re-equilibration). Fractions were collected at 1-min intervals starting from the 10th minute into 10 consecutive tubes (1 → 10 in a cyclic manner), vacuum freeze-dried, and stored at − 80 °C until mass spectrometry analysis.

Liquid chromatography-tandem mass spectrometry (LC–MS/MS): prior to LC–MS/MS analysis, fractionated peptides were reconstituted and spiked with iRT peptides (1:10 ratio) as internal standards. Peptides were analyzed using a high-resolution mass spectrometer (Q Exactive HF-X, Thermo Fisher Scientific) with the following parameters: scan range: 350–1250 m/z; isolation window: 26 m/z; fragmentation: higher-energy collisional dissociation (HCD); resolution: 60,000 (MS1), 15,000 (MS2).

Data processing and bioinformatics analysis: raw spectra were matched against a UniProt human protein database using Spectronaut Pulsar software (Biognosys). Data-independent acquisition (DIA) processing included retention time alignment, peak extraction, and library-based quantification. Missing values were imputed using the k-nearest neighbors (KNN) algorithm, and proteins missing in > 50% of samples or with a coefficient of variation (CV) > 30% were excluded.

Differential protein expression was assessed using one-way ANOVA (p < 0.05). Multivariate analyses, including principal component analysis (PCA) and orthogonal partial least squares-discriminant analysis (OPLS-DA), were performed to evaluate sample clustering and discriminatory features. Functional enrichment analysis was conducted using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. Protein–protein interaction networks were constructed using STRING (v11.5) and visualized in Cytoscape (v3.9.1).

Metabolomics

Fifty microliters of synovial fluid from each sample were mixed with 250 μL of isotope-labeled methanol, vortexed, and incubated at − 20 °C for 20 min. After ultrasonic extraction on ice for 15 min, samples were centrifuged at 13,300 rpm (4 °C, 15 min). The supernatant (200 μL) was collected, vacuum-dried at 30 °C for 2 h, and reconstituted in 0.1 mL of HILIC reconstitution solution. Following sonication, vortexing, and centrifugation, 80 μL of the supernatant was transferred to an LC–MS vial for analysis using an ACQUITY UPLC system coupled to an AB TripleTOF 5600 high-resolution mass spectrometer. For quality control (QC), 20 μL of supernatant from each sample was pooled, aliquoted, and analyzed alongside experimental samples.

Raw data were processed using Progenesis QI (Waters) for peak alignment, annotation, and filtering (retention score ≥ 35, intensity > 1000). After median normalization and missing value imputation (KNN algorithm), PCA assessed data quality. Differential metabolites were identified by t-test (p < 0.05) and OPLS-DA (variable importance in projection (VIP) ≥ 1). Pathway analysis was performed via MetaboAnalyst 4.0 using the KEGG database.

Weighted gene co-expression network analysis (WGCNA)

A total of 376 proteins and 20 samples were analyzed. Proteins with low expression variability (standard deviation ≤ 0.5) were filtered out, leaving 208 proteins and 20 samples for further analysis. A weighted co-expression network model was constructed using a power value of 22, and the remaining 208 proteins were divided into three modules. Data analysis and visualization were performed using the WGCNA package in R, and data visualization was performed using R and Python. The Pearson correlation algorithm calculates the correlation coefficient and p-value between module characteristic proteins and traits. Modules with an absolute correlation coefficient ≥ 0.3 and a p-value < 0.05 were considered significant. For each considerable module, the correlation between module gene expression and trait gene significance (GS) was calculated, and the correlation between module gene expression and Eigengene was analyzed to construct a module-trait correlation analysis.

Integrated omics analysis

Orthogonal partial least squares (O2PLS) analysis was employed to integrate the proteomics data and metabolomics data. O2PLS is a multivariate data integration technique that decomposes the variation in two datasets into joint, dataset-specific, and noise components, enabling the identification of correlations between the datasets while accounting for their intrinsic structure [61].

Prior to analysis, the data were centered around zero and scaled. To determine the optimal number of components, cross-validation was performed using the “crossval_o2m_adjR2()” function, which adjusts for cross-validation in O2PLS analysis. This process yielded values for “n” (number of joint components), “nx” (number of transcriptome-specific components), and “ny” (number of metabolite-specific components). In this study, “n = 3,” “nx = 0,” and “ny = 1” were selected to best fit the omics data to the O2PLS model.

In the O2PLS model, the joint components capture the covariance between the proteins and metabolite data, while the dataset-specific components capture the variation unique to each dataset. The loading values for each variable (gene or metabolite) on the joint components indicate their relative importance in determining the joint variation. Variables with high loading values on the same joint component are strongly correlated. Therefore, by examining the variables with high loading values on the joint components, it is possible to identify groups of proteins and metabolites that are related, potentially reflecting underlying biological processes or pathways.

Malondialdehyde (MDA) assay

Synovial tissue and synovial fluid samples were retrieved from − 80 °C storage. Synovial tissue was homogenized in ice-cold lysis buffer (50 mM Tris–HCl, pH 7.4, 150 mM NaCl, 1% Triton X-100, 1 mM EDTA, and protease inhibitors). After homogenization or lysis, samples were centrifuged at 12,000 × g for 10 min at 4 °C, and the supernatant was collected for analysis. Protein concentration was determined using a BCA protein assay kit (Beyotime, P0009). Reaction setup: blank: 0.1 mL of homogenization buffer, lysis buffer, or PBS. Standard curve: 0.1 mL of serial dilutions of MDA standard (provided in the kit). Samples: 0.1 mL of tissue or synovial fluid supernatant. To each tube, 0.2 mL of MDA detection working solution (Beyotime, S0131S) was added. Samples were mixed thoroughly and heated at 100 °C for 15 min using one of the following methods: Post-incubation processing: samples were cooled to room temperature and centrifuged at 1000 × g for 10 min. Two hundred microliters of supernatant was transferred to a 96-well plate, and absorbance was measured at 532 nm using a microplate reader. A reference wavelength of 450 nm was used for dual-wavelength correction. Quantification: For synovial fluid, MDA concentration (μM) was calculated directly from the standard curve. For tissue samples, MDA content was normalized to total protein (μmol/mg protein). Each sample was analyzed in triplicate (n = 3).

Western blot

Western blotting was performed as described previously [25]. Briefly, synovial samples were lysed and centrifuged at 12,000 g for 25 min at four °C, and the supernatant was collected as the tissue protein solution. Protein concentration was measured using the BCA kit. The primary antibodies used were TNFSF11 (1:1000, ABcloal, A2550); CTSK (1:500, ABcloal, A5871); ADGRE5 (1:2000, ABcloal, A22218); NF-κB (1:1000, ABcloal, A2547); and β-Tubulin (1:1000, ABcloal, AC008). Immunoblots were visualized using BeyoECL Plus (Beyotime, Beijing, China), and protein bands were photographed and stored using the Tanon 2500R gel imaging system (Tanon, Shanghai, China). Band intensity was quantified using ImageJ 1.39 V software (n = 3).

Statistical analysis

All quantitative data are presented as mean ± standard error of the mean (S.E.M.). Statistical analyses were conducted using GraphPad Prism 5.0 (GraphPad Software, San Diego, USA). For comparisons between two groups, a two-tailed unpaired Student’s t-test was employed. For multiple group comparisons, one-way analysis of variance followed by Tukey’s post hoc test was used to adjust for multiple comparisons. A p-value < 0.05 was considered statistically significant.

For proteomic and metabolomic differential expression analyses, all statistical computations were performed in R (version 4.41). Differential expression was conducted using the limma package (or other applicable packages), and raw p-values were adjusted for multiple hypothesis testing using the Benjamini–Hochberg method to control the false discovery rate (FDR). Features with FDR-adjusted p < 0.05 and absolute fold change ≥ 2 were considered significantly altered.

Multivariate statistical analyses, including PCA and O2PLS, were performed using the mixOmics R package. In the O2PLS modeling, variables with VIP scores > 1.0 were considered key contributors to group separation. All data visualization and preprocessing steps were conducted in R using standard packages such as ggplot2 and pheatmap.