The schematic diagram of the experimental processes is shown in Fig. 1. The detailed experimental methods are listed one by one below.
Ethical statement
Animal care and experimental procedures were performed according to the internationally recognized ARRIVE guidelines and guidelines of the Institutional Animal Care and Use Committee (IACUC) of Zhejiang Chinese Medical University. Moreover, the ethics committee of Zhejiang Chinese Medical University approved this study protocol (No. SCXK2016-0010). All experimental procedures were based on the principles of the “three Rs” (Replacement, Reduction and Refinement).
Experimental animals and the establishment of the proteoglycan-induced AS model
Thirty-two locally supplied 3-month-old female Balb/c mice were housed in the AAALAC-accredited Experimental Animal Center at the Zhejiang Chinese Medical University in environmentally controlled conditions: temperature: (25±1) °C, relative humidity: (56±2) % and light (12 h/12 h light-dark cycle). All mice had free access to a standard chow pellet diet and water ad libitum. After a one-week adaptation, all mice were divided randomly into the healthy control group mice (HC group, n=8) and other mice (n=24), which were used for AS modelling. According to our previous published study [22] and the other published literature [28, 29], the AS mice models were induced by i.p. injection with an emulsion of 100μg of cartilage proteoglycans (PG) (Sigma-Aldrich, St. Louis, MO, USA) and 2 mg dimethyldioctadecylammonium (DDA) (Sigma-Aldrich, St. Louis, MO, USA) for 3 times at 21 days intervals. Fourteen weeks after the third i.p. injection (at 20 week), the success of AS modeling was confirmed by measurements of arthritis index (AI) score (each paw had a score over 1) and axial skeleton ankyloses through the animal digital radiography (Aolong, Dandong, China) [22, 28, 30]. Afterward, the confirming AS modeling mice were randomly divided into three groups, including untreated AS mice (AS group, n=8), or AS mice that received moxibustion at acupuncture points (MA group, n=8) or AS mice that received moxibustion at non-acupuncture points (MNA group, n=8) [22, 31].
Moxibustion intervention
Mice in the MA group were fixed on the platform after shaving the hair to expose the skin of moxibustion intervention region. As previously described [22, 32], the moxa cone (height: 120 mm, diameter: 7 mm, “Han Medicine”, Nanyang, China) was lit and held 2 cm above the selecting acupuncture points, namely bilateral BL23 and bilateral ST36 and GV4, for approximately 35 min (including 14 min for each pair of BL23 and ST36 and 7 min for GV4). Previously, according to our previous clinical trial [13], moxibustion at BL23 and ST36 and GV4 acupuncture points have been proven to ameliorate the disease activity and improve the physical function for AS patients. The selecting acupuncture points were located in accordance with the “Veterinary Acupuncture & Moxibustion Atlas” [33]. BL23 is located below the spinous process of the second lumbar vertebra (L2), 1.5cun lateral to the posterior midline. ST36 is located posterolateral to the knee about 5 mm below the fibular head. GV4 is located on the posterior midline, in depression below the spinous process of the second lumber vertebra. The moxibustion intervention was performed for 7 consecutive days in a course with a total of four courses. The non-acupuncture points at the base of the tail were used as controls [34, 35]. In addition, in order to avoid the tail acupuncture points located at distal region of the tail, moxibustion was perform at one-fifth point of tail length from the proximal region of the tail [34, 35]. These non-acupuncture points do not overlap with any other known acupuncture points. In addition to acupuncture point selection, the moxibustion intervention method was performed by the same protocol as the MA group. The mice in the control and model groups were also fixed on a board similar to that in the MA and MNA groups, but did not receive moxibustion stimulation.
Arthritis severity, biological sample collections, histopathology, cytokines and OB related genes
In order to record the effects of moxibustion on arthritis severity in mice paws, the AI score in each mouse was recorded by the same observers, who were blind to the group allocation as previously described [29, 32]. The maximum AI score per paw was 4 and the maximum AI score per mice was 16 [29, 32]. The detailed criteria about the AI scoring system was described in Additional file 1: Table S1. Under anesthesia using isoflurane, blood samples were collected through cardiac puncture. The serum was separated by centrifuging at 3000× g for 15 min and the supernatants were stored at -80 °C. After sacrificing all mice, dissected left ankle joints were fixed with 10% formalin (v/v) for 1 day, decalcified with 5% formic acid (w/v) for 30 days and dehydrated by application of ethanol. After that, dehydrated specimens were processed for paraffin embedding, cut into 5 μm thin slices and stained with hematoxylin-eosin (H&E). Subsequently, these stained slices were imaged with the Virtual Slide Scanner NanoZoomer 2.0-HT (Hamamatsu Photonics, Hamamatsu, Japan). An experienced pathologist, who was blinded to the group allocation, assessed the slides in accordance with cartilage destruction, inflammatory cell infiltration, the narrowing of joint space and synovial hyperplasia as previously described [22, 28, 29]. For measuring cytokines and OB-specific genes’ levels, the serum concentrations of CRP, IFN-γ, IL-17, ALP and OCN were detected using commercial ELISA kits (MEIMIAN, Yancheng, China) following the manufacturer’s instruction.
UHPLC-Q-TOF/MS based metabolomics analysis
Sample Preparation for Metabolomics Study
The freeze serum samples were thawed at 4 ℃ temperature. Next, 100μL serum from each sample was mixed with 400μL precooled solution (methanol: acetonitrile, 1:1, v/v). Then, in order to precipitate proteins, the mixture solution was incubated for one hour at -20 °C after vortex-mixing for 60 s. Subsequently, the solution was centrifugated at 4°C (14,000 × g for 20 min). Lastly, the supernatant was piped into tubes for lyophilization and then mixed with solution containing acetonitrile and H2O (1:1, v/v), which was then ready for UHPLC/Q-TOF-MS analysis.
Chromatography
Metabolomics analysis was performed using an Agilent 1290 Infinity UHPLC system. (Waters Corporation, Milford, USA). The serum chromatographic separation was carried out at 25°C on ACQUITY UPLC BEH Amide (1.7 μm, 2.1 × 100 mm column) and ACQUITY UPLC HSS T3 (1.8 μm, 2.1 × 100 mm column). The flow rate was set at 0.3 mL/min. The optimal mobile phase was composed of eluent A (ddH2O+25 mM ammonium acetate+25 mM ammonia) and eluent B (acetonitrile). The gradient elution programs for the serum were set according to our previous study [22]. The QC samples after every 8 detected samples were used to test the suitability and consistency of the LC-MS system.
Mass spectrometry
The AB Triple TOF 6600 (AB Sciex Corporation, Framingham, MA, USA), operating in both positive and negative ion modes, was applied to electrospray ionization (ESI)-Mass spectrometry (MS) experiments. The key corresponding MS conditions were set as follows: Ion Source Gas1 (Gas1): 40, Ion Source Gas2 (Gas2): 80, Curtain gas (CUR): 30, source temperature: 650℃, Ion Sapary Voltage Floating (ISVF)±5000 V (in ESI+ and ESI- mode); The information-dependent basis (IDA) method based on Analyst TF 1.7 software (AB Sciex Corporation, Framingham, MA, USA) was used for MS acquisition. The mode was set as Declustering potential (DP):±60 V(in ESI+ and ESI- mode), Collision Energy: 35±15 eV. The IDA was set as exclude isotopes within 4 Da, Candidate ions to monitor per cycle:10.
Statistical analysis
The acquired raw UHPLC-Q/TOF-MS data from serum samples were first converted to the common .mzXML format by utilizing ProteoWizard software. Then, these converted datasets were processed using XCMS Online program (http://metlin.scripps.edu) for peak identification and matching, alignment, peak filtration, and translating into the three-dimensional (retention time, mass, intensity of the peaks) data matrices [36]. Subsequently, the resulting dataset of each sample was normalized and transformed by log with Metaboanalyst 4.0 (Montreal, Quebec, Canada) [37]. After that, these resulting data matrices were imported into SIMCA-P 14.0 software package (Umetrics AB, Umeå, Sweden) for a series of multivariate statistical analysis (MVA) after pareto-scaling. Firstly, a partial least squares discriminant analysis (PLS-DA) was performed to overview the metabolic profile differences among four groups. Next, the supervised orthogonal projection to latent structure-discriminant analysis (OPLS-DA) was utilized to distinguish different groups. A permutation test (200 permutations) was used to confirm the validity of the MVA model against over-fitting. Moreover, the quality of the MVA model was evaluated by the R2 (cum) and Q2 (cum) values. The potential metabolic candidates were filtered based on the following: Benjamini-Hochberg’s False Discovery Rate (FDR) correction<0.05 in one-way analysis of variance (one-way ANOVA) followed by Tukey's HSD comparison test using GraphPad Prism 8 software packages (San Diego, California, USA), fold change (FC) >1.33 or <0.77, variable importance for project (VIP) values of the established OPLS-DA model > 1 and the correlation coefficient |r| of the established OPLS-DA model >0.55 [38]. Then, as described in the previous research [39], the enhanced four-dimensional volcano plot including the parameters from the univariate analysis (log2 FC,–log(p-value)) and multivariate analysis (VIP and correlation coefficient (|r|) values of the established OPLS-DA model) were conducted to reflect the potential discriminating variables among the different groups.
Metabolites identification and metabolic pathway analysis
The structure of metabolites was identified by accurate mass number matching (< 25 ppm) and secondary spectrogram matching. The LC/MS spectra was assigned according to the freely accessible online biochemical databases (HMDB, METLIN, KEGG etc.) and in-house developed LC/MS databases.
MetPA (Metabolomics Pathway Analysis) was used for holistic PATHWAY analysis. In the pathway analysis algorithms model, “hypergeometric test” and “relative-betweeness centrality” were selected for performing over-representation analysis and pathway topology analysis, respectively. Moreover, the KEGG category was chosen as the pathway library model. The threshold of the raw P-value was set as 0.05 to filter out the most important and potential relevant metabolic pathways.
DIA-MS based proteomics analysis
Sample Preparation and Fractionation for DDA library Generation
Firstly, in order to effectively expand and improve the dynamic range, 150 µl of pooled serum from each group was immunodepleted for removal of top 14 high abundant proteins using a commercial depletion kit (MARS™, Agilent Technologies) coupled to a High-Performance Liquid Chromatography (HPLC, Shimadzu LC-10AT). High abundant proteins depletion was performed according to the manufacturer’s protocols and buffer exchange was employed with 25 mM ammonium bicarbonate (ABC) in ddH2O using ultrafiltration devices with a 10 kD molecular weight cut off (Microcon, Millipore-Merck, Germany). After collecting the resulting protein solution, the total protein concentrations were quantified with a Bicinchoninic acid protein assay (Thermo, IL, USA) according to the manufacturer’s instructions.
Secondly, the serum protein digestion steps were carried out following the Filter-Aided Sample Preparation approach [40]. In short, protein concentrates (200 μg) were lysed directly using 30 μl of STD lysis buffer (4% [w/v] SDS, 100 mM Dithiothreitol in 0.15 M Tris pH 8.0). Afterward, excess reagents in the resulting lysates of serum samples were depleted using UA solution (8 M Urea in 0.15 M Tris pH 8.0) by repeated ultrafiltration (30 kDa cut-off, Microcon, Millipore-Merck, Germany). After that, 100 μl of UA buffer containing 0.05 M iodoacetamide (IAA) was mixed with the filters, and the sample mixtures were incubated at room temperature for 30 min in darkness. Subsequently, filter units were washed three times with 100 μl of UA buffer and then twice with 100 μl of 25 mM Ammonium Bicarbonate (Sigma, MO, USA). Finally, each filter unit was sequentially digested with 2 μg trypsin (Promega, Madison, WI, USA) in 40 μl of 100 mM Ammonium Bicarbonate buffer and the mixture solution was incubated overnight at 37 °C. All resulting tryptic peptides were then collected via centrifugation. Peptide content was estimated by UV light spectral density at 280 nm.
Thirdly, pooled tryptic peptides were next fractionated using a commercial peptide fractionation kit (Pierce™ High pH Reversed-Phase Peptide Fractionation Kit, Thermo Fisher Scientific). Each fraction was concentrated by vacuum centrifugation (Benchmark Scientific, NJ, USA) and further acidified with 10µl of 0.1% (vol/vol) formic acid. Afterwards, all resulting peptides were acidified with 40µl of 0.1% (v/v) formic acid prior to desalting with C18 Cartridges (Empore™ SPE Cartridges C18 (standard density), bed I.D. 7 mm, volume 3 ml) according to the manufacturer’s protocol (Thermo Fisher Scientific, MA, USA) [41].
Fourthly, the iRT-Kits (Biognosys, Schlieren, Switzerland) were added into the pooled digest serum samples to normalize the relative retention time.
DDA Mass Spectrometry Assay
All fractions were used for generating an extensive serum-specific DDA library. As previously described, DDA Mass Spectrometry Assay was carried out with the help of a Thermo Scientific Q-Exactive HF mass spectrometer connected to an Easy-nLC™ 1200 chromatography equipment (Thermo Scientific, MA, USA) [41]. Firstly, two micrograms of peptides were re-dissolved in 2% (v/v) ACN and 0.1% (v/v) FA, and then pooled tryptic peptide samples were injected into an EASY-SprayTMC18 Trap column (Thermo Scientific, P/N 164946, 3um, 75um*2cm). Subsequently, a 120-minute segmented gradient of solvent [80% (v/v) ACN and 0.1% (v/v) FA] at a flow rate of 250 nl/min on an EASY-SprayTMC18 LC Analytical Column (Thermo Scientific, ES803, 2um,75um*50cm) was used for peptide separation. Afterwards, eluting peptides were detected by mass spectrometry analysis in an ESI+ ion model. The parameters for one MS survey scan experiment were set as follows: (1) The mass range: m/z 300–1,650 Da; (2) High resolution, R: 60,000 (FWHM) at m/z 200; (3) Dynamic exclusion: 30s; (4) AGC (Automatic gain control) target for MS: 3×106; (5) Maximum IT: 25 ms. Each full MS SIM scan followed 20 dd-MS2 scans, and the parameters for tandem MS/MS scans experiments were set as follows: (1) High resolution, R: 15, 000 (FWHM) at m/z 200; (2) AGC target for MS2: 5×104; (3) Maximum IT: 25 ms; (4) Fragmented by normalized collision energy (NCE): 27 eV.
Mass Spectrometry Assay for DIA
Each sample peptide was analyzed by LC-MS/MS performing in the DIA mode with the Q-Exactive HF mass spectrometer (Thermo Scientific, MA, USA) by Shanghai APT proteomic institute (Shanghai, China). In the DIA analysis, each DIA cycle comprised of one full MS–SIM scan with m/z 15 Da per isolation window to cover the scan range of m/z 350–1,650 Da [42]. The experimental set-up conditions for the DIA analysis were almost similar to the DDA analysis. The MS1 scans were set as: (1) High resolution, R: 60,000 (FWHM) at m/z 200; (2) AGC target for MS: 3×106; (5) Maximum IT: 50 ms and profile mode. MS2 scans were set as: (1) High resolution, R: 30, 000 (FWHM) at m/z 200; (2) AGC target for MS2: 3×106; (3) Maximum IT: auto; (4) Fragmented by NCE: 30 eV. Moreover, equal aliquot of each sample from all groups were pooled as QC samples were injected with DIA mode at the beginning of the LC-MS/MS analysis. The QC samples after every 4 detected samples were used to examine the suitability and consistency of the LC-MS system.
Mass spectrometry data analysis
For protein quantification, DDA library data files were processed with MaxQuant analysis software (version 1.5.3.17), which performed a tandem mass spectrometry (MS/MS) spectra search against the combined databases of FASTA sequences and iRT standard peptide sequences. Searching parameters were set as follows: (1) Enzyme specificity: full tryptic//P digestion; (2) Control modification: carbamidomethylation of cysteine (C); (3) Dynamic modification: N-acetylation (N-term) of proteins and oxidation of Methionine (M); (4) Tryptic peptides: no more than 2 missed cleavages; (5) The satisfy FDR: ≤ 1% at peptide level. For generating a reference DDA spectral library, the original raw files and DDA searching results of each sample were imported into the Spectronaut Pulsar XTM (version 12.0.20491.4) software (Biognosys, AG, Switzerland). The building DDA spectral library was next normalized to iRT standard peptides.
Subsequently, Spectronaut Pulsar XTM (version 12.0.20491.4) software (Biognosys, AG, Switzerland) was applied to match the original DIA-MS data against the above constructed DDA spectral library. The parameters were set as follows: (1) RT prediction type: Dynamic iRT; (2) Interference correction on MS2 level; (3) Cross run normalization. All results were quantified and filtered based on adjusted Q value < 0.01 criterion.
Protein identification and bioinformatics methods
The DEPs were filtered by Benjamini-Hochberg’s FDR correction<0.05 in one-way ANOVA followed by Tukey's HSD comparison test using GraphPad Prism 8 software packages (San Diego, California, USA) as well as FC >1.33 or <0.77 [32].
In order to determine the functional classification and functional enrichment of DEPs altered in the serum of AS mice by moxibustion interventions, Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.8, an online functional annotation tool, was adopted to perform a Gene Ontology (GO) enrichment analysis. Biological Networks Gene Ontology (BiNGO), a plug-in in a cytoscape visualization platform software (version 3.7.2, free available from http://www.cytoscape.org/), was employed to further annotate functional Go terms into three categories, including biological functions, molecular function and subcellular locations [43]. Canonical pathway enrichment analysis of DEPs was carried out using the Kyoto Encyclopedia of Genes and Genomes (KEGG) (https://www.kegg.jp/) automatic annotation server and REACTOME (https://reactome.org/) databases. FDR values less than 0.05 were considered to indicate statistically significant enriched GO terms or associated canonical pathways. The Protein-Protein interaction (PPI) network was established by STRING online software.
Validation of DEPs
Four of the DEPs identified by the DIA-MS based proteomic analysis (LCAT, TF, HSPA8 and MIF) were validated using commercial enzyme-linked immunosorbent assay (ELISA) kits (MEIMIAN, Yancheng, China) according to the manufacturer’s instructions. The OD value was read at 450 nm with a PerkinElmer EnSpireTM multilabel reader (PerkinElmer Life and Analytical Science, Turku, Finland).
Integrated analysis of metabolomics and proteomics
A multi-omic analysis was implemented for integrating the results of the single proteomic and metabolomic. According to the data of proteomics and metabolomics platform, MetaboAnalyst 4.0 was carried out to acquire the integrated canonical pathway. The significantly enriched integrated pathways were filtered after examining a threshold of FDR <0.05. The correlation network diagram established by OmicBean software (a web-based software application for interacting multi-omic data and available from http://www.omicsbean.cn/) can display moxibustion-intervened differentially expressed metabolites, DEPs, and significantly perturbed signaling pathways through visual interaction. The protein–metabolite interaction network was generated with a score of ≥ 2, number of nodes>60 and p-value < 0.05 [44]. Finally, the hierarchical cluster analysis of identified proteins and metabolites was employed by MEV software (version 4.9.0).
Afterward, the Pearson correlation matrix analysis approach through Correlation Calculator (version 1.0.1) was performed to explore the correlations between the potential biomarkers including validated proteins, significant metabolites and biochemical indicators [45]. A value of p<0.05 was considered statistically significant.