Multi-omics Analysis at Serum Proteomic and Metabolomic Levels Reveals Mechanisms of Moxibustion for a Mouse Model of Ankylosing Spondylitis


 Background: Moxibustion is widely used in ameliorating symptoms of ankylosing spondylitis (AS). The aim of this study was to identify the potential molecular profiles of moxibustion on relieving AS mice through incorporating ultra-high-performance liquid chromatography quadrupole-time of flight mass spectrometry (UHPLC-Q-TOF/MS) based metabolomics approach and data independent acquisition-mass spectrometry (DIA-MS) based proteomic. Methods: In this research, mice with proteoglycan-induced spondylitis (PGISp) were intervened with moxibustion for four weeks at specific acupuncture points or at non-acupuncture points. Arthritis severity, histopathological examinations, cytokines and osteoblast (OB) related genes were assessed. Protein profiling and metabolite profiling of serum of mice were examined by UHPLC-Q-TOF/MS based metabolomics technology and DIA-MS based proteomic approach. A multi-omic analysis was implemented for integrating the results of the single proteomic and metabolomic. The levels of some novel proteins in significantly enriched integrated pathways were validated by the enzyme-linked immunosorbent assay (ELISA) method. Results: Our findings clearly indicate that acupuncture points’ moxibustion can significantly decrease the arthritis index (AI) score, improve the histopathological examination and reduce the serum levels of C-reactive protein (CRP), Interferon-γ (IFN-γ), Interleukin-17 (IL-17), Alkaline phosphatase (ALP) and Osteocalcin (OCN) for AS mice. A total of 21 differentially expressed metabolites and 21 differentially expressed proteins (DEPs) were identified as potential moxibustion-intervened biomarkers for AS mice. These molecules were mainly involved in six integrative canonical pathways: mineral absorption, lipid metabolism, purine metabolism, glycolysis/gluconeogenesis, glycine, serine and threonine metabolism and phenylalanine metabolism pathways. Furthermore, The ELISA analysis of four key proteins was consistent with the results of the proteomic analysis. Conclusions: Combing UHPLC-Q-TOF/MS based metabolomics approach and DIA-MS based proteomic, as a promising tool, can advance our understanding of potential mechanisms of action of moxibustion for AS from a holistic perspective.

experienced pathologist, who was blinded to the group allocation, assessed the slides in accordance with cartilage destruction, in ammatory cell in ltration, the narrowing of joint space and synovial hyperplasia as previously described [22,28,29]. For measuring cytokines and OB-speci c 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 H 2 O (1:1, v/v), which was then ready for UHPLC/Q-TOF-MS analysis.

Chromatography
Metabolomics analysis was performed using an Agilent 1290 In nity 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 ow rate was set at 0.3 mL/min. The optimal mobile phase was composed of eluent A (ddH 2 O+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.

Statistical analysis
The acquired raw UHPLC-Q/TOF-MS data from serum samples were rst 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 identi cation and matching, alignment, peak ltration, 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 pro le 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 con rm the validity of the MVA model against over-tting. Moreover, the quality of the MVA model was evaluated by the R 2 (cum) and Q 2 (cum) values. The potential metabolic candidates were ltered 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 coe cient |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 (log 2 FC,-log(p-value)) and multivariate analysis (VIP and correlation coe cient (|r|) values of the established OPLS-DA model) were conducted to re ect the potential discriminating variables among the different groups.

Metabolites identi cation and metabolic pathway analysis
The structure of metabolites was identi ed 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 lter 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 ddH 2 O using ultra ltration devices with a 10 kD molecular weight cut off (Microcon, Millipore-Merck, Germany). After collecting the resulting protein solution, the total protein concentrations were quanti ed 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 ultra ltration (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 lters, and the sample mixtures were incubated at room temperature for 30 min in darkness. Subsequently, lter 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 lter 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 Scienti c). Each fraction was concentrated by vacuum centrifugation (Benchmark Scienti c, NJ, USA) and further acidi ed with 10µl of 0.1% (vol/vol) formic acid. Afterwards, all resulting peptides were acidi ed 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 Scienti c, 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-speci c DDA library. As previously described, DDA Mass Spectrometry Assay was carried out with the help of a Thermo Scienti c Q-Exactive HF mass spectrometer connected to an Easy-nLC™ 1200 chromatography equipment (Thermo Scienti c, 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 Scienti c, 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 ow rate of 250 nl/min on an EASY-Spray TM C18 LC Analytical Column (Thermo Scienti c, 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)  Subsequently, Spectronaut Pulsar X TM (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 quanti ed and ltered based on adjusted Q value < 0.01 criterion.

Protein identi cation and bioinformatics methods
The DEPs were ltered 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 classi cation 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 signi cant 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 identi ed 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 signi cantly enriched integrated pathways were ltered 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 signi cantly 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 identi ed 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, signi cant metabolites and biochemical indicators [45]. A value of p<0.05 was considered statistically signi cant.

Results
AI score, the histological examinations, and expression of in ammatory cytokines and OB related genes In line with the literature [28,30], fourteen weeks after the third i.p. injection (at week 20), mice in the AS, MA and MNA group showed signi cant increasing levels of AI score when compared with the HC group (p<0.01). Before performing the moxibustion (at week 20), there were no statistically signi cant differences among the AS, MA and MNA group in terms of the levels of AI score (p>0.05). However, compared with the AS group, the AI score in the MA group was signi cantly decreased after 4-week moxibustion therapy (p<0.05). ( Fig. 2A).
As shown in Additional le 1: Figure S1, histopathological results of the left ankle joint from AS model group (Figure S1 B) showed cartilage enlargement and erosion, in ammatory cell in ltration and narrowing the joint space compared with those in the control group. Moreover, compared with the AS group, both MA group and MNA group can remarkably inhibit the narrowing of joint space in AS mice ( Figure S1C and Figure S1D). However, compared with the MNA group, there is more signi cant improvement in synovial hyperplasia, in ammatory cell in ltration and cartilage defection of AS mice with moxibustion at speci c acupuncture points (Figure S1 C). It suggests that moxibustion intervention at speci c acupuncture points had a relatively better effect on ameliorating the arthritic characteristics in AS mice.
As shown in the Fig. 2 B-F, the serum levels of CRP, IFN-γ, IL-17, ALP and OCN in the AS group were signi cantly increased as compared with the HC group (p<0.01). However, compared with the AS and MNA group, moxibustion on the GV4, BL23 and ST36 acupuncture points signi cantly decreased the serum levels of CRP, IFN-γ, IL-17, ALP and OCN (p<0.01). Compared with the AS group, moxibustion on non-acupuncture points had favorable but not statistically signi cant effects in reducing the serum levels of CRP, IFN-γ, IL-17, ALP and OCN (p>0.05).

UHPLC/Q-TOF-MS method validation
In order to eliminate system technical errors, the relatively clustering of QC samples was calculated to examine the repeatability and robustness of the UHPLC/Q-TOF-MS metabolomics method. As shown in Additional le 1: Figure S2, all of the four QC samples in both ESI+ and ESI-ion modes were less than 2 SD, indicating that the data from the UHPLC/Q-TOF-MS metabolomics were reliable [46].

Pattern recognition analysis
The representative UHPLC/Q-TOF-MS total ion chromatogram (TIC) pro les of serum samples in ESI+ and ESI-mode are shown in Additional le 1: Figure S3.
Although the TICs were quite different from the ESI+ to ESI-mode, visual inspection of TICs in the same ESI+ or ESI-mode of HC, AS, MA and MNA groups showed less clear differences. Thus, MVA, including PLS-DA and two-group of OPLS-DA, was used to further discover any potential endogenous metabolites related to AS and moxibustion treatment.
PLS-DA was rst utilized to overview the metabolic pro le differences of serum samples among different groups. PLS-DA of serum dataset showed a clear group separation among HC, AS, MA and MNA groups in both ESI+ and ESI-modes (Fig. 3). Moreover, according to Fig. 3, variations of serum metabolic pro ling in the MA group were much closer to the HC group than others; whereas scattered points from the MNA group were much closer to the AS group than the MA group. Therefore, metabolic perturbations caused by the AS modeling had the tendency to be reversed moving towards to the normalized state following moxibustion at speci c GV4, BL23 and ST36 acupuncture points, but little in uenced by non-acupuncture points moxibustion intervention.
Identi cation of potential serum biomarkers and the changing trends among different groups Herein, pairwise comparisons using the OPLS-DA model were further performed to elucidate the potential biomarkers that were responsible for the pathological process of AS and possible anti-AS mechanism of moxibustion intervention in serum samples. After pareto scaling, OPLS-DA score plot showed good separations between HC vs. AS ( Fig. 4A and 4B) and AS vs. MA ( Fig. 4C and 4D) in both ESI+ and ESI-mode. The cumulative R 2 Y and Q 2 of OPLS-DA model in ESI+ and ESI-modes suggested that the models were good to prediction and reliability [47] (HC vs. AS: ESI+: R 2 Y=0.998 cum, Q 2 =0.870 cum; ESI-: Additional le 1: Figure S4B and S4D). Then, the enhanced four-dimensional volcano plot revealed a number of metabolic features that were responsible for the differences among different groups. In this study, 61 discriminatory endogenous metabolites were ultimately screened out which showed signi cant alterations between HC and AS groups according to the pre-de ned criteria. Notably, MA intervention remarkably regulated 21 of 61 metabolites (Benjamini-Hochberg's FDR correction<0.05, AS vs. MA in one-way ANOVA followed by Tukey's HSD comparison test) ( Fig. 4E and Table 1). Histograms of all 21 differential metabolites in both ESI+ and ESI-mode were shown in Additional le 1: Figure S5, which included 2 elevated metabolites and 19 decreased metabolites. Additionally, 9 metabolites (phosphorylcholine, glucose 6-phosphate, L-Phenylalanine, allantoin, serine, ethanoic acid, D-Mannose, L-Malic acid and citrate) were signi cantly altered in both MA and MNA groups compared with the AS group, most of which were associated with the TCA cycle and energy metabolism (Additional le 1: Figure S5 and Table 1). Finally, the clustering heat-map showed the trends of the serum metabolic concentration difference among different groups (Fig. 4F).

Metabolic pathway analysis
As shown in Fig. 5, 21 endogenous metabolites listed in Table 1 primarily involved in Glyoxylate and dicarboxylate metabolism, Aminoacyl-tRNA biosynthesis, Phenylalanine, tyrosine and tryptophan biosynthesis, Phenylalanine metabolism, Glycine, serine and threonine metabolism, Glycerophospholipid metabolism, Citrate cycle (TCA cycle), Purine metabolism and Pyruvate metabolism.

DIA-MS based Proteomics results
Expressions of moxibustion-associated proteins in the serum of AS mice DIA-MS based proteomic platform was used to comprehensively pro le the moxibustion-trigged protein expression in the serum of AS mice. To minimize biological variation, two or three serum samples were pooled together to acquire three biological replicates for each group. Ultimately, DIA-MS based proteomic analysis, identi ed a total of 7,614 peptides, corresponding to 486 proteins across pooled serum samples from the HC, AS, MA and MNA groups. Then, 439 of 486 proteins were successfully quanti ed for further statistical analysis. Following the pre-de ned criteria, 34 DEPs were regarded as potential proteomic pro les of PG-induced AS with FDR <0.05 (HC vs. AS in one-way ANOVA followed by Tukey's HSD comparison test), and FC>1.33 or <0.77. Importantly, 21 of 34 proteins in the serum of AS mice were signi cantly reversed after moxibustion intervention at GV4, BL23 and ST36 acupuncture points (FDR <0.05, AS vs. MA in one-way ANOVA followed by Tukey's HSD comparison test), which included 6 up-regulated proteins and 15 down-regulated proteins ( Table 2). The hierarchical clustering of these 21 candidate proteins revealed that the HC, AS, MA and MNA groups were considerably variable from each other (Fig. 6). While, only 4 of 21 proteins in the serum of AS mice (PSMA3, APOB, XDH and PGAM1) were signi cantly altered in both MA and MNA groups compared with the AS group, most of which were involved in the proteasome pathway and energy metabolism (Table 2).
Next, GO annotation was further performed to elaborate the sophisticated biological functions, molecular functions and subcellular localization of 21 target proteins of moxibustion intervention applied to the serum of AS mice. Totally, 631 biological functions terms, 111 molecular function terms, and 69 subcellular localization terms were signi cantly enriched after screening with a threshold of FDR<0.05. The top 10 most enriched Go terms for each category are shown in Fig. 7. Notably, most of the enriched biological function terms were concentrated in the lipoprotein metabolic process and in ammatory response (Fig. 7A).
The molecular function terms affected by moxibustion intervention primarily involved in threonine peptidase activity, lipid and cholesterol metabolic process, xanthine oxidase activity and protein binding (Fig. 7C). In terms of the subcellular localization category, extracellular space, lipoprotein particle and proteasome complex comprised the most enriched terms (Fig. 7E). Then, these top 10 statistically signi cant GO terms in each category (biological functions, molecular functions and subcellular localization) and their relationship with adjacent functional terms in Biological Networks were visually by BiNGO, a plug-in Cytoscape software. (Fig. 7B, 7D and 7F).
Furtherly, in order to gain insights on the pathway alterations associated with the moxibustion intervention, the KEGG enrichment analysis was performed with the aid of Omicbean software. The results showed that multiple canonical pathways were altered, and most of them were associated with the metabolic pathways (10%), cholesterol metabolism (8%) and proteasome (6%) (Fig. 8A). Moreover, the bubble plot further revealed that the statistically signi cant canonical pathways that include 21 DEPs contributing to the effect of moxibustion were cholesterol metabolism, proteasome, ferroptosis, legionellosis, caffeine metabolism, sulfur metabolism, ubiquinone and other terpenoid-quinone biosynthesis, phenylalanine metabolism, protein processing in endoplasmic reticulum and vitamin digestion and absorption (Fig. 8B). Finally, the Protein-Protein interaction (PPI) network of 21 altered proteins affected by moxibustion intervention is depicted in Additional le 1: Figure S6.
Integrated pathway and network analysis and validation of DEPs via ELISA For better explorations on the moxibustion-targeted molecular information, we integrated metabolomics and proteomics data, to obtain shared canonical pathways and interaction networks between these signi cantly altered small molecular metabolites and proteins with MetaboAnalyst 4.0 and OmicBean software. Overall, 27 signi cantly integrated pathways were identi ed with p<0.05, whereas 7 canonical pathways comprised both proteins and metabolites ( Fig. 9A and Table S2). Interestingly, most of altered canonical pathways were related to the metabolic pathways, which are somewhat in line with the results from the single proteomics or metabolomics analysis. Our previous study used a 2 DE combined with MALDI-TOF-MS/MS technique to explore potential biomarkers of moxibustion a CIA rat model [32]. The results suggested that the metabolic pathway was the most important biological pathway involved in the regulation of moxibustion, which was also somewhat consistent with the present study.
Notably, according to the Fig. 9B, LCAT, the only detected protein in the glycerophospholipid metabolism pathway, was also involved in another enriched pathway: cholesterol metabolism pathway. Similarity, glycine, serine and threonine metabolism pathway and glycolysis or gluconeogenesis pathway shared the common protein: PGAM1. Moreover, XDH, which has a direct link to the several metabolites (uric acid, xanthine and hypoxanthine) that included in the purine metabolism pathway, may also play a vital role in the predicted protein-metabolite network. Moreover, TF and HSPA8 have the greatest number of close connections and interactions with the downstream proteins of the mineral absorption pathway. Specially, HSPA8 has the closest relationship with the PSMB2, PSMA6, and PSMA3, which are the crucial proteins involving in the proteasome pathway. MIF, the only detected protein in the phenylalanine metabolism pathway, can independently predict radiographic progression in the spine of AS. Thus, the above important proteins may play a domain role in this predicted protein-metabolite network. However, only LCAT, TF, HSPA8 and MIF are meridian speci c proteins (Table 2). Thus, the above four key proteins in the predicted network were chosen to be further validated by ELISA. The ELISA results revealed that four DEPs were signi cantly up-regulated (LCAT and HSPA8) or downregulated (TF and MIF) in the serum of the MA group compared with the AS model group (Fig. 10), consistent with DIA-MS based proteomic analysis.

Correlation between potential biomarkers and biochemical indicators
The results of correlation analysis between potential biomarkers and biochemical indicators in serum were shown in Additional le 1: Figure S7. At the protein level, TF was positively and highly correlated with IFN-γ, IL-17 and AI score and MIF had a remarkable positive correlation with CRP. However, LCAT was negatively correlated with IFN-γ in our research. Moreover, TF and MIF showed obvious positive correlation with all OB-speci c genes including ALP and OCN.
At the metabolic level, there was a signi cant positive correlation between taurine and cytokines such as CRP and IL-17. However, tyrosine and D-Mannose were negatively related to the in ammatory measurements. Moreover, L-phenylalanine was closely related to the OB-speci c gene ALP.

Discussion
Similar to our previous studies [22,32], the ndings of this study clearly indicate that acupuncture points' moxibustion can signi cantly improve the in ammatory status of AS mice, which includes decreasing the AI score and reducing the serum levels of CRP, IFN-γ, IL-17 in AS mice. Moreover, moxibustion at speci c GV4, BL23 and ST36 acupuncture points also improves histopathological examinations and decreases the serum levels of OB-speci c genes, including ALP and OCN for AS mice. According to results of metabolomic and proteomic studies, 21 differentially expressed metabolites and 21 DEPs were identi ed as potential moxibustion-intervened biomarkers for AS mice. Furthermore, the predicted protein-metabolite network showed that six altered canonical pathways, including mineral absorption, lipid metabolism, purine metabolism, glycolysis/gluconeogenesis, glycine, serine and threonine metabolism and phenylalanine metabolism pathways, were closely related to the anti-AS effect of moxibustion. The details will be further discussed accordingly.
Previously, several acute-phase reactants, including serotransferrin (TF), are over-expressed in the PBMCs of AS patients compared with healthy controls, and anti-TNF-α agents signi cantly decreased the levels of TF in AS patients [48,49]. In the present study, the serum level of TF increases in AS modeling mice compared with HC mice, and it decreases markedly after 4-weeks' moxibustion intervention, which is in line with the previous clinical studies. Thus, moxibustion may also have a direct effect on decreasing the acute-phase protein levels in AS. Meanwhile, TF is an iron-regulator protein which can mediate the uptake of iron in many cells. The previous research had revealed that there was a negative correlation between the levels of iron ion and the serum TF in Sprague-Dawley rats [50]. In this study, the signi cantly decreasing expression of TF that is regulated by moxibustion intervention may contribute to the increasing the concentrations of iron, which is consistent with results from the anti-TNF therapy in the previous clinical report [49]. Using vitro murine calvarial MC3T3-E1 OB cells culture system, Yamasaki et al. [51] had revealed that high concentrations of iron could prohibit the proliferation, differentiation and deposition of calcium in MC3T3-E1 OB cells. Notably, AS is speci cally characterized by OB activity and bone resorption increasing [52,53]. Thus, moxibustion may suppress the OB activity, calci cation and mineralization for AS mice by reducing the serum levels of TF, which is further veri ed by the down-regulation of the OB speci c gene (ALP and OCN). (Fig. 10 and Additional le 1: Figure S7) In addition to TF, other acute-phase proteins in the mineral absorption pathway such as Hp, Cp and ORM1, were also up-regulated in AS modeling mice. Hp plays a signi cant role in modulating immune and in ammatory responses on various immune cells (neutrophils, macrophages, lymphocytes, etc.).
Previously, Li et al. had reported an up-regulation of Hp levels in patients with AS, but the serum levels of Hp are not strongly correlated with the disease activity of AS [54,55]. Cp is a kind of metal proteins, involving in excretion of copper, copper-dependent oxidation, and iron transporting [56]. A great number of studies had already reported that the serum level of Cp was over-expressed in patients with rheumatic conditions (including AS) compared with healthy individuals, which was consistent with the results of our current study [57,58]. ORM1 is a transport protein which is involved in regulating the activity of the immune system during the acute-phase reaction. Additionally, combined serum levels of ORM1, Apcs, and CRP were shown to have a high value for diagnosing AS [59]. However, moxibustion intervention can apparently reverse the over-expression levels of these three acute-phase proteins, suggesting that moxibustion plays a critical role in decreasing the acute phase process of AS.
Hemopexin (HPX), primarily a heme binding protein, is frequently secreted during in ammatory responses. Using iTRAQ-based proteomic approach, Cai et al. [28] had demonstrated that compared with the healthy volunteers, the HPX protein was found to be over-expressed in PBMCs of AS patients, which was in line with the results of our current study. To explore this issue, HPX was known as a 'toxic heme' transporter and scavenger in the blood, thus the up-regulating of HPX in AS could protect cells against oxidative stress [60,61]. In this research, the serum level of HPX is further elevated after moxibustion at GV4, BL23 and ST36 acupuncture points, suggesting that moxibustion is able to promote the host defense and inhibit oxidative stress for AS. Also, another previously published proteomic study from our research team had supported a signi cant increase in the HPX concentration after moxibustion in the collagen-induced arthritis (CIA) rat [62].
HSPA8 is a type of heat shock proteins (HSPs). The upregulation expression of HSPs may exert a protective effect on organisations by stabilizing the protein structure, inhibiting the release of proin ammatory cytokines (TNF-α, IL-1β, IL-6, etc.) and modulating the host immune response, especially during the chronic stress condition [63]. Of note, HSPs may also serve as crucial roles involved in the mechanism of action of moxibustion [64]. Previously, Yi et al [65] and Chang et al [66] proved that moxibustion at ST36 and ST21 acupuncture points could protect gastric mucosal cells against stress injury by increasing the levels of HSPs. Moreover, our previous proteomic study [32] also had proven that moxibustion at ST36 and BL23 acupuncture points could elevate the serum levels of HSPs, which played a key role of anti-arthritis in the CIA rat. Consistent with these previous reports, in this study, the serum concentration of HSPA8 has increased markedly following moxibustion intervention at GV4, BL23 and ST36 acupuncture points, and AI score reduction ( Fig. 2A) and ankle pathological improvement (Additional le 1: Figure S1) also have been observed in the MA group than in the AS model group, implying the protecting effect of moxibustion on musculoskeletal damage in AS mice.
Based on our studies, three main downstream enzymes (PSMB2, PSMA6 and PSMA3) of HSPA8, were found to be over-expressed in the serum of AS mice.
Notably, these proteins are involved in the proteasome pathway that we have found in the proteomic KEGG enrichment analysis (Fig. 8B). These results, somewhat, consist with that of Write et al [67] who reported that proteins involved in the proteasome pathway were up-regulated in the monocytes of AS patients. To explore this issue, proteasome is a proteolytic enzyme complex existing in the cytoplasm and nucleus, which is responsible for degrading most proteins in the cell, and plays a key role in the expression of MHC class I molecules (including B27 molecules) and antigen presentation [68]. However, the activation of proteasome could degrade MHC class I molecules produce antigen arthritic peptides that can bind to their own B27 molecules, and further present to CD8 cytotoxic T cells, thereby initiating the pathogenesis of AS. In the present study, moxibustion at proper acupuncture points could prohibit the activation of proteasome pathway via downregulating the expression of PSMB2, PSMA6 and PSMA3 and played a role of anti-AS [69].
Citric acid is a vital intermediate of the TCA cycle and the VCP is primarily involved in the ATP metabolic process. Both of them play a prominent role in the metabolism of energy. It had been well reported that metabolites and proteins involved in the TCA cycle were signi cantly raised to enhance energy production, so as to compensate for insu cient energy sources caused by glycolysis under AS in ammation condition [70]. Similarity, in this study, the levels of citric acid and VCP increased in the PG-induced AS model mice group compared with the healthy control group, suggesting a boosted energy supply to compensate for the shortage of body energy. Previously, our research team had demonstrated that moxibustion could correct abnormal energy metabolism in CIA through down-regulating the plasma levels of VCP [32]. Moreover, the urine level of citric acid has also been proven to down-regulate signi cantly following moxibustion intervention in AS mice [22]. Consistent with these ndings from our previous research, herein, by means of acupuncture point's moxibustion intervention, decreased serum levels of citric acid and VCP are observed in this study, which indicate that moxibustion at speci c acupuncture points can contribute to reducing the symptoms of AS by recovering the impairment of energy metabolism.

Lipid metabolism pathway
In this research, the molecular alterations in the serum, including lipoproteins (APOA2, APOB and APOC1), LCAT, and LysoPCs (Sn-Glycerol 3phosphoethanolamine; Glycerophosphocholine; Phosphatidylcholine), were mainly associated with the lipid metabolism pathway. Interestingly, as shown in Fig. 9B, the cholesterol metabolism and glycerophospholipid metabolism which were two important parts of the lipid metabolism pathway shared one common lipoprotein, LCAT. Previously, a great number of reports had revealed that AS patients with high disease activity was closely related to a deteriorated lipid pro le, which suggested that lipid metabolism pathway might be involved in the pathogenesis of AS [71,72].
Compared with the HC group, proteins associated with cholesterol metabolism including APOA, APOB, APOC1 and LCAT were differently regulated in the AS model group. APOA is a high-density lipoprotein (HDL) apolipoprotein, which can transfer cholesterol from tissues to the liver for decomposition. Importantly, Lecithin cholesterol acyltransferase (LCAT) is the key enzyme in this metabolic process, which can elevate the expression of APOA and further promote the clearance of cholesterol [73]. However, in this study, in ammatory conditions induced by AS decreased capacity of HDL and LCAT to e ciently transfer excess cholesterol in the ligament tissue to the liver, leading to the ectopic fat deposition in the ligaments of AS model mice [74]. Then, ectopic fat deposition could further promote the production of in ammatory cytokines, which resulted in cartilage dysfunction and degeneration in patients with AS. Our previous LC-MS based tissue metabolomics research [22] had shown that moxibustion intervention could correct the abnormal ectopic fat deposition through decreasing the tissue metabolic levels of cholesterol in AS model mice. The results of this multi-omic study further support our previous metabolomic study, suggesting that moxibustion at speci c acupuncture points can ameliorate the ectopic fat deposition of AS by up-regulating the expression levels of APOA and HDL proteins related to cholesterol synthesis and transport.
In contrast, the other two lipoproteins, including APOB and APOC, signi cantly increased in the AS model mice. These results were in agreement with our previous clinical proteomic study, which showed that APOB and APOC were over-expressed in the serum of AS patients compared with healthy controls [75].
APOB and APOC were the main protein components of low-density lipoproteins (LDL) and very low density lipoproteins (VLDL), respectively, and they were important biomarkers for predicting the cardiovascular (CV) risks [76]. It has been reported that the CV event is the leading cause of mortality in AS, which is partly due to the dyslipidemia in AS patients, including high levels of APOB and APOC [77]. Previously, Huang et al [78] reported that plasma levels of APOB and APOC markedly decreased in the hyperlipidemia rabbit model and increased following herbal-partitioned moxibustion. In line with the previous research, in this study, moxibustion may also well have a clinically relevant effect on decreasing the CV risks of AS by down-regulating the serum levels of APOB and APOC.
LysoPCs, the main active components of oxidized low-density lipoprotein (ox-LDL) in the glycerophospholipid metabolism pathway, have been con rmed to be generated by lecithin after being hydrolyzed by LCAT. It was reported that high concentrations of LysoPCs could active neutrophils and macrophage and cause oxidative stress damage in the cell membrane [79]. Previously, moxibustion had exerted a protective effect on cell membrane of rats with ethanol-induced gastric mucosal lesions through down-regulating the serum levels of LysoPCs [31]. Similarity, the levels of three LysoPCs compounds are reversed by moxibustion at GV4, BL23 and ST36, suggesting that moxibustion at speci c acupuncture points treatment has the ability to ameliorate the cellular injury in AS.

Purine metabolism
It was reported that the purine metabolism was notably involved in the pathogenesis of adult and pediatric spinal spondyloarthritis, and overproduction of proteins and metabolites in purine metabolism could aggravate joint and cartilage tissue injury for AS patients [80,81]. Consistent with those previous studies, the up-regulation of XDH and its downstream metabolic products inosine, hypoxanthine, xanthine, xanthosine and uracil in AS model mice in the current study also suggested that AS had critically in uenced purine metabolism. XDH is a key rate-limiting enzyme in purine degradation, which contributes to transferring hypoxanthine to xanthine. Inosine and uracil are downstream metabolites of nucleotides through the purine metabolism pathway. Moreover, inosine can also be served as novel small-molecule fecal biomarkers in AS and showed a close relationship with the in ammatory indicator ESR [82]. Previously, methotrexate (MTX), the second line disease modifying antirheumatic drug prescribed in AS patients, exerted anti-in ammatory effects for AS via inhibiting purine metabolism [83]. Moreover, Du et al [84] performed the LC-Q/TOF-MS metabolomics approach to discover the effect of heat-reinforcing acupuncture at ST36 on urinary metabolic ngerprinting of arthritis rabbit with cold syndrome. They also found that compared with the model group, the purine metabolism in the heat-reinforcing acupuncture group was signi cantly decreased. Results of this study exhibited that the molecular levels of XDH, inosine, hypoxanthine, xanthine, xanthosine and uracil in AS mice were markedly down-regulated by acupuncture points' moxibustion therapy, suggesting a similarity of anti-AS mechanisms to MTX and acupuncture.

Glycolysis/gluconeogenesis
Glycolysis/gluconeogenesis, another important metabolic pathway that supplies ATP for organisms, plays a key role in maintaining the energy homeostasis.
PGAM1 is an important metabolic enzyme that participates in the glycolysis/gluconeogenesis pathway by interconverting 3-and 2-phosphoglycerate with 2,3bisphosphoglycerate. It has been reported that the expression levels of the PGAM1 in PBMCs of AS patients are higher than those of healthy volunteers, implying that glycolysis is activated under systemic in ammatory conditions, which further results in a hypoxic microenvironment that induces the chondrocyte apoptosis [85,86]. However, in this study, moxibustion can ameliorate the hypoxic microenvironment and reduce chondrocyte apoptosis in AS mice via down-regulating the serum levels of PGAM1 that involve in the glycolysis/gluconeogenesis pathway, which is in agreement with our previous studies [22,32]. Moreover, PGAM1 changes were also closely related to the improvement of pathological results of cartilage tissue (Additional le 1: Figure S1) and the decrease of AI score in the MA group ( Fig. 2A).
Glucose 6-phosphate, D-Mannose and ethanoic acid, the downstream metabolites of PGAM1, were all markedly up-regulated in the AS model mice, which was in line with our proteomic study and further supported the enhancing activity of glycolysis and dysfunction of energy metabolism in AS. Our previous urinary metabolomic study [22] had revealed that metabolites that were involved in glycolysis/gluconeogenesis were all signi cantly down-regulated following moxibustion intervention. In this study, moxibustion may correct the impairment of energy metabolism through decreasing the serum levels of Glucose 6phosphate, D-Mannose and ethanoic acid, which was somewhat in line with the previous research.

Glycine, serine and threonine metabolism
The altered proteins and metabolites involved in the glycine serine and threonine metabolism has been reported in PBMCs isolated from the plasma of AS patients [85] and the fecal extracts of AS modeling rats [87]. These past ndings are consistent with our current results.
Notably, PGAM1 is not only involved in glycolysis/gluconeogenesis but also simultaneously participates in glycine, serine and threonine metabolism. PGAM1 has a lot of subfamily protein members, including Bisphosphoglycerate Mutase (BPGM), which can control glycolytic intermediate levels and further modulate the serine synthesis. Thus, the up-regulation of serine concentration caused by the increasing of PGAM1 levels in glycolysis/gluconeogenesis pathway can result in the misfolding of the HLA-B27 protein process, which plays an important role in the pathogenesis of AS [88]. However, the protein level of PGAM1 together with the metabolic concentration of serine in the moxibustion at acupuncture points group are lower than that in the AS model group, which suggests that moxibustion can simultaneously improve the energy homeostasis and ammonia acid metabolism.
Glycine, as a decomposition product by serine, has been revealed to possess anti-in ammatory properties [89]. Thus, lower levels of glycine in serum samples of AS model mice and the depletion is reversed by 4-week moxibustion at speci c acupuncture points, indicating that AS modeling may cause a weaken antiin ammatory properties and acupuncture points' moxibustion therapy has a positive effect in improving this trend.
Taurine, as a sulfur-containing AA, plays an important role in directly scavenging reactive oxygen species and protecting joint tissue from oxidative stress in progressive in ammatory diseases [90]. In this research, the serum level of taurine signi cantly down-regulates in AS model mice suggesting that the dysfunction of taurine and hypotaurine metabolism and decreasing in antioxidant capacity exist during the pathological process of AS. Recently, Zhang et al [91] proved that moxibustion at ST36 and BL23 acupuncture points could up-regulate the serum level of taurine in CIA rats. Similarity, the decreasing AI scores and increasing concentration of taurine following moxibustion at speci c acupuncture points in this study seemed to be a serum biomarker of the AS model mice in local defense reactions against oxidative stress injury.

Phenylalanine metabolism
During the development of AS, the dysregulation of tight junctions between intestinal epithelial cells increased gut permeability then led to subsequent gut leaky. Pathogenic microbiota from gut leaky of the intestinal epithelium could further result in the dysregulation of osteoproliferation in bone marrow through activating innate immunity system and IL-23/Th17 axis of AS [92]. Interestingly, the Macrophage migration inhibitory factor (MIF) is also mainly involved in the innate immune response of pathogenic microbiota caused by gut leakage. Some speci c gut pathogens, including the Helicobacter pylori, have been proven to trigger the release of MIF [93]. Subsequently, the production of MIF from the gut may drive the new bone formation in patients with AS via activating and mineralizing OBs, upregulating genes involved in osteogenesis, and triggering the stabilization of β-catenin, an important mediator involved in the OBspeci c Wnt/β-catenin pathway [70]. Thus, it is not surprising to detect the high concentration of MIF in patients with AS [94]. Moreover, the serum level of MIF could also independently predict radiographic progression in the spine of AS patients [95]. However, moxibustion may directly decrease the osteogenesis of AS through downregulating the serum levels of MIF, which is in line with the reducing trends of ALP and OCN after moxibustion intervention ( Fig. 10 and Additional le 1: Figure S7).
Increased levels of L-Phenylalanine, tyrosine and L-histidine in the AS model group could re ect early changes in colonic protein fermentation and the disturbances of intestinal homeostasis, which may further result in abnormal metabolism of intestinal ora and phenylalanine [87]. The previous study had proved that herbal-partitioned moxibustion could protect the integrity of intestinal barrier in rats with ulcerative colitis by reversing microbiota pro les and rebalancing the intestinal microbial environment [96]. Moreover, our previous research [97] had found that moxibustion at ST36 and BL23 acupuncture points could attenuate in ammatory symptoms of CIA model rats to some extent via decreasing the urine level of L-phenylalanine and tyrosine. In line with previous works, serum levels of L-Phenylalanine, tyrosine and L-histidine are down-regulated after moxibustion intervention at GV4, BL23 and ST36 acupuncture points in AS mice, which indicate that moxibustion can ameliorate the dysbiosis of intestinal ora and the perturbed phenylalanine metabolism pathway in AS.

Limitation
There are several potential limits of this study. Firstly, we only select a single time point to explore the possible molecular mechanism of moxibustion for AS mice. However, the time-relative therapeutic mechanism of moxibustion on AS is still scarce. Secondly, some lipoproteins and LysoPCs were found to be involved in the anti-AS effects of moxibustion based on the DIA-MS based proteomic and UHPLC-Q-TOF/MS based metabolomics. In the future, targeted lipidomics technique is encouraged to further elucidate the possible lipid metabolism mechanism of moxibustion for treating AS. Thirdly, moxibustion may exert anti-AS effect through alleviating the dysbiosis of intestinal ora and the perturbed phenylalanine metabolism pathway. However, considering the complexity of gut microbiota, the intestinal ora mechanism of moxibustion in the treatment of AS can be also further studied in-depth by combining with 16S rDNA gene sequencing and metagenomic sequencing technologies.

Conclusions
In conclusion, moxibustion achieved its anti-AS effect probably by modulating mineral absorption, lipid metabolism, purine metabolism, glycolysis/gluconeogenesis, glycine, serine and threonine metabolism and phenylalanine metabolism pathways at system level. Combing DIA-MS based proteomic and UHPLC-Q-TOF/MS based metabolomics approach, as a promising tool, can advance our understanding of potential mechanisms of action of moxibustion for AS from a holistic perspective.     Graph of the top 10 signi cant pathways of biological function category (A), molecular function category (C), and subcellular localization category (E). The blue color frame indicates the statistical differences of the GO terms via adjust P-value using Benjamini-Hochberg method. A cut-off red line represents the threshold (adjust p value=0.05). After log transformation, p<0.05 will have the value of -log (p-value)>1.33. The pink color frame represents the number of proteins. Enrichment analysis of the biological function category (B), molecular function category (D), and subcellular localization category (F); each node means one different biological function, molecular function, or subcellular localization and more signi cant ones are lled in with a deeper color.