Proteomic profiling analysis of postmenopausal osteoporosis and osteopenia identifies potential proteins associated with low bone mineral density

Study populations and design

For this study, we recruited a total of 54 postmenopausal females, including 18 female patients with primary PMOP, 18 with postmenopausal ON, and 18 with normal bone mass as NCs were recruited from March 2018 to February 2019 in Honghui Hospital, Xi’an Jiaotong University (People’s Republic of China). Patients with secondary OP, including those with Cushing’s syndrome, hyperparathyroidism, hyperthyroidism, steroid abuse, and other malignancies such as chronic liver disease, heart disease and kidney disease were excluded. Drug factors including glucocorticoids, immunosuppressants, heparin, anticonvulsants, anticancer drugs, aluminum-containing antacids, thyroid hormones, etc. were also excluded. In addition, BMD measurements were not available in the following situations: (1) those orally administered with drugs within 2–6 days of measurement that affect image development; (2) when radioisotope inspections were carried out within 3 days; (3) those who cannot lie on the examination bed, or cannot adhere to the 5-min test; (4) those with a severe deformity of the spine or metal implants on the spine. Finally, biochemical detection of serum 25-hydroxyvitamin D (vitamin D), alkaline phosphatase (ALP), procollagen type 1 N-peptide (P1NP) and C-terminal cross-linking telopeptide of type 1 collagen (CTX) concentrations as well as BMD measurements were determined for all participants. Each participant signed an informed consent form, and the project was approved by the Institutional Review Board, Honghui Hospital, Xi’an Jiaotong University (Project No: 2018-22). We collected participants’ information, including their age, age at menopause, height (m) and body mass (kg). Participants’ BMI was calculated as mass/height². The experimental workflow is presented in Fig. 1.

BMD measurement

Dual-energy X-ray absorptiometry was used to detect the BMD of the lumbar spine (L_1–4) and hip (Discovery Wi, Hologic, Marlborough, MA, USA) in all participants. The scan images were obtained from the detector and are shown in Fig. 2. We collected information pertaining to patients’ BMD values (g/cm²) and T scores. The T-score reference ranges were calculated using data from healthy young Asian women provided by the bone densitometry equipment manufacturer. Participants were categorized as having OP, ON, or normal BMD according to the World Health Organization’s T-score classification: participants with T-scores ≤–2.5 at any site were diagnosed as having OP, –2.5 < T-scores < –1.0 had ON and T-scores ≥ –1.0 were indicative of normal bone mass. The PRM approach has been used for BMD-related DAPs validation.

Blood sample collection

We collected the fasting venous blood samples from each participant, and the blood samples were centrifuged at 4 °C at 3,000 rpm for 15 min. After the centrifugal, the supernatant from each sample were collected and stored at −80 °C until measurement. In each group, we randomly mixed six plasma samples into one sample for the TMT proteome experiment and PRM validation.

Biochemical detection

Serum 25-hydroxyvitamin D (Vitamin D), alkaline phosphatase (ALP), procollagen type 1 N-peptide (P1NP), and C-terminal cross-linking telopeptide of type 1 collagen (CTX) were detected using the Roche electrochemiluminescence system (Cobas e 601 analyzer, Roche Diagnostics GmbH, Wetzlar, Germany).

Protein extraction and enzymolysis

The plasma pools were depleted of the proteins by using the Multiple Affinity Removal LC Column-Human 14 (Cat. 5188-6560, Agilent, CA, USA) according to the protocol from the manufacturer. The desalination and concentration of the low-abundance components were performed by using a 10 kDa ultrafiltration tube. One volume of SDT (4% SDS, 100 mM Tris/HCl pH 7.6, 0.1 M dithiothreitol) buffer was added, and then boiled for 15 min, then centrifuged for 20 min with 14,000×g. The BCA Protein Assay Kit (Bio-Rad, Hercules, CA, USA) was used for protein quantification, and all the supernatant samples were stored at −80 °C for further measurement.

SDS–PAGE separation

Overall, proteins (20 µg) were mixed by the 5× loading buffer for each sample, respectively, and boiled for 5 min. Proteins from all the samples were separated by a 12.5% sodium dodecyl sulfate (SDS)-polyacrylamide gel electrophoresis (PAGE) gel and the Coomassie Blue R-250 staining system was used for binds visualization.

TMT labeling

A total of 200 μg of proteins for each sample were incorporated into 30 μl SDT buffer (4% SDS, 100 mM DTT, 150 mM Tris-HCl pH 8.0). The detergent, DTT and other low-molecular-weight components were removed using UA buffer (8 M Urea, 150 mM Tris-HCl pH 8.0) by repeated ultrafiltration (Microcon units, 10 kD). Then 100 μl iodoacetamide (100 mM IAA in UA buffer) was added to block reduced cysteine residues and the samples were incubated for 30 min in darkness. The filters were washed with 100 μl UA buffer three times and then 100 μl 100 mM TEAB buffer twice. Finally, the protein suspensions were digested with 4 μg trypsin in 40 μl TEAB buffer overnight at 37 °C, and the resulting peptides were collected as a filtrate. The peptide content was estimated by UV light spectral density at 280 nM using an extinctions coefficient of 1.1 of 0.1% (g/l) solution that was calculated on the basis of the frequency of tryptophan and tyrosine in vertebrate proteins.

Following this and according to the manufacturer’s instructions (Thermo Fisher Scientific, Waltham, MA, USA), a peptide mixture of each sample (100 μg) was labeled with the TMT reagent.

Peptide fractionation with high pH reversed phase

According to the manufacturer’s instructions, increasing the acetonitrile step-gradient elution, the TMT-labeled digest samples into 10 fractions fractionated with a Pierce high pH reversed-phase fractionation kit (Thermo Fisher Scientific, Waltham, MA, USA).

Mass spectrometry

Each fraction was injected in order to analyze the nano-LC-MS/MS. A reverse-phase trap column which is connected to a C18 reversed-phase analytical column in buffer A (0.1% formic acid) and separated with a linear gradient of buffer B (84% acetonitrile and 0.1% formic acid) loaded the peptide mixture and the flow rate was 300 nL/min and it was controlled by IntelliFlow technology.

The LC-MS/MS analysis was completed with a Q Exactive mass spectrometer and it was coupled to Easy nLC for 60 min. The mass spectrometer was performed in positive ion mode. Using a data-dependent top10 method can obtain the MS data, and the most abundant precursor ions were selected from the scan of 300–1,800 m/z. An automatic gain control (AGC) target was 3e6, and the maximum inject time was 10 ms. Dynamic exclusion duration was 40 s. Survey scans were obtained with the resolution of 70,000 at m/z 200 and the resolution for the HCD spectra was 17,500 with m/z 200, 35,000 at m/z 200, and the width of isolation was 2 m/z. The normalized collision energy was 30 and the under fill ratio specified the minimum percentage of the target value was defined as 0.1%. The instrument was performed by using peptide recognition mode enabled. The LC–MS/MS spectra were searched by using a MASCOT engine which was embedded into Proteome Discoverer 1.4. The MS proteomics data have been deposited to the ProteomeXchange Consortium via the PRIDE (Perez-Riverol et al., 2019) partner repository with the dataset identifier PXD017804.

PRM analysis

We further performed LC-PRM analysis (Peterson et al., 2012) for several selected proteins from the same samples, to test the protein abundance obtained by previous TMT analysis. We prepared the peptides and then analyzed by using TMT analysis according to the protocol, and TOMHAQ has used for validation (Shanghai Applied Protein Technology Co., Ltd., Shanghai, China) (Zhang et al., 2019). The raw data were analyzed with the bioinformatic tool Skyline (MacLean et al., 2010) (MacCoss Lab, University of Washington, Seattle, WA, USA), where the signal intensities of the peptide sequences for the significantly altered proteins were quantified relative to each sample, and normalized to the standard reference (Liu & Lv, 2019).

Statistical and bioinformatics analysis

Basic information including participants’ age, age at menopause, duration of menopause, BMI, biochemical detection, BMD and T score were displayed as the mean ± standard deviation (SD). A P-value < 0.05 with two-tailed Student’s t-test was considered statistically different between two groups. A correlation analysis between basic information, protein expression, and BMD values was performed using R 3.6.0 tools and GraphPad Prism 8.0.1 statistical software (GraphPad, La Jolla, CA, USA). Proteins with a fold change (FC) > 1.2 as well as a statistical P-value < 0.05 between two groups were selected as DAPs. The unsupervised principal component analysis (PCA) and partial least squares discriminant analysis (PLS-DA) model of proteome data were used in this study. PCA and PLS-DA analysis of proteome data, and the bi-clustering and volcano plots of DAPs were performed with R 3.6.0 tools. We searched Blast2GO (version 3.3.5) during the gene ontology mapping and annotation procedure (Gotz et al., 2008). GO enrichment on three modules including biological process (GO-BP), molecular function (GO-MF), and cellular component (GO-CC) was applied based on a Fisher’ exact test, with P-values < 0.05 considered as statistically significant functional categories.

Characteristics of the study population and correlation analysis

Our study recruited a total of 54 postmenopausal females, including 18 females with normal bone mass (normal controls (NC)), 18 females with ON and 18 with OP. We collected and analyzed patients’ information, including their age, age at menopause (M_age), body mass index (BMI), biochemical detection results and BMD measurement results (BMD value and T score data from the DXA–BMD detection). Figure 2 displays the DXA scan images obtained during the BMD measurement using the detector. As for BMD and T scores, when comparing the ON vs. NC, OP vs. NC, and OP vs. ON groups, all showed statistically significant differences, which is consistent with the diagnostic criteria for OP and ON. There were no statistically significant differences between any of the two groups in terms of age, age at menopause, and serum concentration of Vitamin, ALP, P1NP and CTX, while the participants in OP group showed a significant lower BMI compared with that in NC and ON groups (Table 1).

We further performed Pearson’s correlation analysis between the different indicators to explore BMD-related factors. Overall, the BMD values were negatively correlated with age (Fig. 3A). Meanwhile, a negative correlation was also observed between the menopause duration and BMD values (Fig. 3C), while a positive correlation was found between BMD and age at menopause (Fig. 3B), and between BMD and BMI (Fig. 3D). Besides, we observed a positive correlation between the concentration of Vitamin D and BMD values, although not statistically significant (Fig. S1A), and the concentration of ALP, P1NP, CTX was not significantly related with the BMD values (Figs. S1B–S1D).

DAP analysis among OP, ON, and NC with proteomic data

Before the proteomics experiment, we randomly mixed six samples into one sample in each group; that is, three mixed plasma samples in each group were prepared for further proteome analysis. High-throughput, label-free, liquid chromatography–tandem mass spectrometry (LC–MS/MS) experiments identified a total of 1,092 proteins in the plasma samples in our study. The unsupervised principal component analysis (PCA) and partial least squares discriminant analysis (PLS-DA) model of proteome data revealed evidence of the separation between low BMD (OP and ON) and the NC, especially when discriminating between OP and NC samples, while it showed less-than-ideal effects on the separation between the OP and ON samples (Figs. 4A and 4B). We selected the proteins with a FC > 1.2 (i.e., where upregulation was >1.2 times or downregulation <0.83 times), and P < 0.05 indicated a statistically significant difference in the DAPs between two compared groups. When compared with the NC group, there were 17 upregulated and 20 downregulated DAPs in the ON group (Figs. 4C and 4D), and 32 upregulated and 123 downregulated DAPs in the OP group (Figs. 4C and 4E). Compared with the ON group, there were 11 upregulated and 5 downregulated DAPs in the OP group (Figs. 4C and 4F). If we set the FC > 1.5 (i.e., when the upregulation was >1.5 times or downregulation <0.67 times), P < 0.5 as the criteria for DAPs, there were 8, 127 and 3 DAPs in the ON vs. NC, OP vs. NC and OP vs. ON groups, respectively (Fig. 4C). The details of the protein quantification results and differential analysis between each group comparison are presented in Supplemental Files 1 (ON vs. NC), Supplemental Files 2 (OP vs. NC) and Supplemental Files 3 (OP vs. ON). The peptide lists of all the identified proteins in the study are presented in Supplemental Files 4.

These results may suggest that the proteome of OP samples were different for the NC group, whereas the ON group (the middle group) showed greater similarity with the OP group, as there was a partly overlapping PLS-DA model between the OP and ON groups and the number of DAPs was less in the OP vs. ON group when compared with the ON vs. NC group. Thus, in subsequent analyses, we placed greater focus on the difference between low BMD (OP and ON) and NC. Functional analysis indicated that the DAPs in the ON vs. NC groups mainly demonstrated catalytic, binding, and molecular function regulator activity; the DAPs appeared to be involved in metabolic, cellular, stimulus response, biological regulation and immune system processes (Fig. 5A). Similarly, the DAPs in the OP vs. NC groups mostly exhibited binding, molecular function regulator, transporter, and molecular transducer activity, and appeared to be involved in cellular, stimulus response, biological regulation, metabolic and biological regulation processes (Fig. 5B).

Identification of candidate proteins correlated with low BMD values in PMOP

In our study, most of the DAPs existed in the OP vs. NC groups, followed by the ON vs. NC and OP vs. ON groups. We focused on the DAPs in the low BMD group, which was found in both the OP vs. NC and ON vs. NC comparisons. When compared with the normal BMD samples, seven proteins were differentially expressed in the OP and ON groups, including Adipocyte plasma membrane-associated protein (H0Y512), Protein disulfideisomerase A5 (Q14554), cDNA FLJ51711 (B4DE30), Glucosidase (A0A024R592), Lysozyme C (P61626), Uncharacterized protein DKFZp666N164 (Q658S4) and Complement C4-A (P0C0L4) (Table 2). Figure S2 showed the MS/MS spectra for the seven DAPs in the low BMD groups. Bi-clustering analysis of these seven DAPs indicated a clear discrimination between normal BMD and low BMD samples (Fig. 6A). The detailed expression of the seven proteins in each group were showed in Figs. 6B–6H. The expression of protein protein Lysozyme C was positively related with BMD values (Fig. S3A), while the expression of proteins proteins Glucosidase, cDNA FLJ51711, and Protein disulfideisomerase A5 was negatively related with BMD values (Figs. S3B–S3D). The correlation between the expression of Complement C4-A, Uncharacterized protein DKFZp666N164, Adipocyte plasma membrane-associated protein, and BMD values was not highly significant (Figs. S3E–S3G). The reported numbers of proteins were supported with more than one peptide.

Verification of the candidate protein biomarkers associated with low BMD with PRM

In this study, Adipocyte plasma membrane-associated protein, Protein disulfideisomerase A5, cDNA FLJ51711, Glucosidase, Lysozyme C, Uncharacterized protein DKFZp666N164 and Complement C4-Awere differentially expressed in OP and ON group compared with normal sample. However the correlation between the expression of Complement C4-A, Uncharacterized protein DKFZp666N164, Adipocyte plasma membrane-associated protein, and BMD values was not highly significant. Therefore we further selected the four BMD-related DAPs including Lysozyme C, Glucosidase, cDNA FLJ51711, and Protein disulfideisomerase A5 for validation using the PRM approach. Among them, the abundance of Lysozyme C was also decreased in the low BMD groups, while Glucosidase and Protein disulfideisomerase A5 expression were increased as BMD decreased (Fig. 7). These results from the PRM testing suggested that the candidate proteins demonstrated similar trends as the TMT results, which confirmed the credibility of the proteomics data and verified the DAPs between the NC and low BMD samples.