Metabolomic analysis of biochemical changes in the tissue and urine of proteoglycan-induced spondylitis in mice after treatment with moxibustion.

BACKGROUND: Moxibustion is widely used in East Asian countries to manage the symptom of rheumatic diseases. The aim of this study was to identify potential metabolic profiles of moxibustion on relieving ankylosing spondylitis (AS) mice through UHPLC-Q-TOF/MS metabolomic study.
METHODS: Thirty-two female Balb/c mice were randomized into healthy control (HC), AS model, moxibustion at acupuncture points (MA) in AS model, and moxibustion at non-acupuncture points (MNA) AS model groups. Moxibustion was administered daily at GV4, bilateral BL23 and bilateral ST36 acupuncture points for four weeks in the MA group. The overall health status, the thickness of hind paws and the tissue concentrations of IL-1β, PGE2, IL-6 and TNF-α were assessed. The UHPLC-Q-TOF/MS was used to explore the perturbations of endogenous metabolites in tissue and urine of AS model mice intervened by moxibustion.
RESULTS: Compared with the AS group, the overall health status was significantly improved after 4-week moxibustion intervention (p < 0.05). The results also showed that MA significantly reduced the levels of paw thickness and decreased the levels of four cytokines in the tissue (p < 0.01). Thirty-seven endogenous metabolites identified by the OPLS-DA were considered to be contributing to therapeutic effects of moxibustion. Moreover, metabolic pathway analysis further revealed that the identified metabolites were mainly involved in TCA cycle, Lipid metabolism, Amino Acid metabolism, Intestinal flora metabolism and Purine metabolism.
CONCLUSIONS: UHPLC-Q-TOF/MS based metabolomics approach, as a novel and powerful tool, can help us to gain the insights into potential mechanisms of action of moxibustion for AS.

Introduction

Ankylosing spondylitis (AS) is a chronic rheumatic disease of unknown exact causes and typical clinical manifestations of AS can include inflammatory back pain, stiffness, deformity and ankylosis of the spine and etc. Recent epidemiological survey shows that about 0.2–0.3% of the Chinese population is suffering from AS, of which young people aged between 20 and 30 years appear to be most frequently affected.

ASAS/EULAR international guidelines suggested nonsteroidal anti-inflammatory drugs (NSAIDs), and subsequently systematic biological agents therapies. However, these recommended NSAIDs medicines are frequently associated with cardiovascular risk increasing and adverse gastrointestinal and renal effects. Although several dramatic biological agents have shown their benefits on reducing the activity for moderate-to-severely AS patients but its’ cost is high. Therefore, Chinese AS patients have turned their attention to seek traditional Chinese medicine (TCM). Moxibustion is one of the most popular therapies for treating rheumatic diseases (including AS) in China.

Several clinical trials and systematic reviews suggested that moxibustion as an adjuvant therapy in combination with NSAIDs or biological agents may exert positive effects on AS patients.6, 7, 8 Potential mechanisms have been suggested that inhibition of non-specific endochondral ossification of spine, regulating immune function, anti-inflammation, anti-oxidant, anti-muscle degenerative and cartilage-protective effect.

Metabolomics strategies are fully complied with the holistic and multi-target characteristics of TCM and have been successfully applied to explore the biochemical mechanisms of moxibustion therapy in the treatment of irritable bowel syndrome (IBS), gastric mucosal lesions (GML), and chronic atrophic gastritis (CAG). Nevertheless, there have been no studies specifically focusing on the metabolic alterations corresponding to moxibustion interventions for AS. Therefore, this study was investigated the mechanism of moxibustion for treating AS with a quantitative metabolomics approach based on ultra high performance liquid chromatography quadrupole-time of flight mass spectrometry (UHPLC-Q-TOF/MS) combined with a multivariate statistical analysis.

Methods

Mice and the creation of the proteoglycan-induced AS model

All animal experiments were performed according to the U.S. National Institutes of Health ‘Guide for the Care and Use of Laboratory Animals’ and guidelines of the Institutional Animal Care and Use Committee of Zhejiang Chinese Medical University. Moreover, this study protocol was also approved by the ethical review of experimental animal welfare in Zhejiang Chinese Medical University (No. SCXK2016-0010). All susceptible 3-month-old female Balb/c mice purchased from Cavens Animal Lab (Hangzhou, China) were housed in the AAALAC-accredited SPF grade Animal Research Center at the Zhejiang Chinese Medical University (Hangzhou, China) in an environment with controlled temperature (25 °C), humidity (56%) and light (12 h light-dark cycle) with uninterrupted access to a standard chow pellet diet and water ad libitum. The animals were acclimatized for 1 week before use. Afterward, eight mice were randomly chosen as the healthy control mice (HC group, n = 8) while the other mice were reserved for AS modeling. As previously well described , 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 swelling of the peripheral feet via the digital vernier caliper (Gongxing, Shanghai, China) and axial skeleton ankyloses through the animal digital radiography (Aolong, Dandong, China). Then, a computer-generated randomized code was used to assign the confirming AS modeling mice to three different 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).

Moxibustion protocols

After AS modeling, mice from both MA and MNA groups were subjected to moxibustion intervention using moxa sticks with a diameter of 7 mm and length of 120 mm (Han Medicine Co Ltd, Nanyang, China). In the MA group, moxibustion was performed at GV4, bilateral BL23 and bilateral ST36 acupuncture points, located in accordance with the “Veterinary Acupuncture & Moxibustion Atlas”. The detailed information regarding the location of GV4, BL23 and ST36 acupuncture points as well as an acupuncture point diagram in mice are summarized in Supplementary A. These identifying acupuncture points were based on TCM theory and found effective in ameliorating symptoms of AS in our previous clinical trial. Prior to moxibustion stimulation, mice in the MA group were fixed in a prone position on a board after shaving the hair to expose the skin and sterilizing the skin surface at GV4, bilateral BL23 and bilateral ST36 acupuncture points. As previously described , the moxa cone was lit for about 5 seconds and then placed on choosing acupuncture points with the fire head 2 cm from the skin surface. Moxibustion at each acupuncture point lasted about 7 min and the total time of moxibustion intervention for each mouse lasted approximately 35 min (including 14 min for each pair of BL23 and ST36 and 7 min for GV4). Each mouse was given daily moxibustion intervention by TCM nursing practitioner for 4 weeks. The tail was used as a non-acupuncture point moxibustion control site (Supplementary A). In addition to acupuncture point selection, the moxibustion intervention method was performed by the same protocol as the MA group. The mice in the HC and AS groups were also fixed on a board similar to that in the MA and MNA groups, but did not receive moxibustion stimulation.

The overall health status measurements

The overall health status of the mice was monitored according to the body weight, physical activity, and behavioral observation scale. (1) The body weights of all mice were recorded in grams by using the electronic weighing balance (Kaifeng, Jinhua, China). (2) The physical activity was measured by gait score, which graded from 0 to 3 as follows: 0: normal gait; mice run and walk normally, 1: slight disability; mice run and walk with difficulty, 2: moderate disability; mice walk with difficulty due to intermittent loading of inflamed joint, and 3: severe disability; three- legged gait. The gait score was performed in a blinded manner and always carried out by the same experimenter, to minimize variability. (3) Behavioral observation scale and its three domains (fur, mental state, and behavior) as a specific instrument to measure the overall health status of the mice (Supplementary B). The score for each domain was provided on a scale ranging from 0 to 3, thus, the total score of this scale was ranged from 0 to 9, with lower scores indicating a better overall health status.

Physical parameters, biological sample collections and cytokines

The maximal thickness of hind paws in each mouse was measured by the digital vernier caliper (Gongxing, Shanghai, China). The samples of 24 hour urine in each mouse were collected in ice-cooled Eppendorf tubes containing 1% NaN3 after 4-week moxibustion intervention. After that, animals were sacrificed under isoflurane anesthesia and ligament tissue samples of the spine were excised and snap-frozen in liquid nitrogen. All above urine and tissue samples were stored at -80 °C immediately for further analysis. For measuring cytokine levels, the tissue concentrations of IL-1β, PGE2, IL-6 and TNF-α were detected using commercial ELISA kits (R&D, USA) following the manufacturer's instruction.

Sample Preparation for Metabolomics Study

The freeze-dried tissue samples were thawed at ambient temperature. Next, 1000 μL of precooled solution (methanol: acetonitrile: H2O, 2:2:1, v/v/v, as an internal standard, concentration: 1 mg/ml) and 50 mg of tissues were mixed. After vortex-mixing for 30 s, the mixture solution was homogenized for 4 minutes with the help of an automatic rapid grinding machine (Jingxin Technology, Shanghai, China). Then, in order to precipitate proteins, the homogenized solution was followed by sonication for 10 min in an ice-water bath, incubation for one hour at -20 °C and centrifugation at 4 °C (12,000 g for 15 min). Subsequently, 750 μL of the supernatant was transferred into a new micro-centrifuge tube that was evaporated to dryness via a vacuum concentrator. Lastly, the remaining dry extract was reconstituted in 100 μL of mixture solution containing acetonitrile and H2O (1:1, v/v), vigorous shaked for 30 seconds, sonicated for 10 min in an ice-water bath and centrifuged at 12,000 g for 15 min. Then, a total of 60 μL of the supernatant was collected into a 2 ML LC/MS glass vial, which was then ready for UHPLC/Q-TOF-MS analysis. Moreover, 10 μl aliquot of each supernatant from all groups of the samples were pooled as the quality control (QC) sample for further UHPLC-Q/TOF-MS analysis. The QC samples after every 7 detected samples were used to test the suitability and consistency of the LC-MS system. For urine samples, the urine samples were thawed at ambient temperature, and 200 μL of each urine sample were chosen and handled by the same method as tissue samples.

Metabolomics Analysis Based on UHPLC-Q/TOF-MS

Chromatography

Metabolomics analysis was performed using an ACQUITY 1290 UHPLC system from Waters Corporation (Waters Corporation, Milford, USA). Tissue and urine chromatographic separation was performed at 40 °C on an Acquity UHPLC BEHC Amide C18 column (1.7 μm *2.1mm*100 mm). The optimal mobile phase was composed of eluent A and eluent B, which were 25 mM ammonium acetate and 25 mM aqueous ammonia in water (A) and 0.1% formic acid in acetonitrile (B), respectively. The gradient elution programs for tissue and urine are shown in the below: (0→0.5 min, eluent A 5% and eluent B 95%; 0.5→7 min, eluent A 5%→35% and eluent B 95%→65%; 7→8 min, eluent A 35%→60% and eluent B 65%→40%; 8→9 min, eluent A 60% and eluent B 40%; 9→9.1 min, eluent A 60%→5% and eluent B 40%→95%; and 9.1→12 min, eluent A 5% and eluent B 95%).

Mass spectrometry

The AB 5600 Triple TOF Mass spectrometer (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: Normalized collision energy: 30 eV; Atomizing gas: Nitrogen; Nebulizer gas: 60 Psi; Ion source gas: 60 Psi; Curtain gas (CUR): 35 Psi; temperature (TEM): 650°C; Ionspray voltage floating (ISVF): 5000 v in positive ion modes (ESI + ) and -4000 v in negative ion modes (ESI-). Information-dependent basis (IDA) method based on Analyst TF 1.7 software (AB Sciex Corporation, Framingham, MA, USA) was used for MS acquisition. In each acquisition cycle, the most intense ions (mass ranges greater than 100 Da) were chosen for the fragmentation (15 MS events per 50 ms).

Data Handling and statistical analysis

The acquired raw UHPLC-Q/TOF-MS data in instrument-specific format (.d) 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. After that, these resulting data matrixes 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 utilized to overview the metabolic profiles differences of tissue and urine samples in the HC, AS, MA and MNA group. Next, the supervised orthogonal projection to latent structure-discriminant analysis (OPLS-DA) was carried out to verify the PLS-DA model, and further maximize distinguish and separate metabolic alterations among groups. The fitness and predictability of the MVA model were controlled and explained by the R2Y (cum) and Q2 (cum) values, respectively. The potential metabolic candidates were filtered based on the following: adjusted p-value (False Discovery Rate, FDR) <0.05 in one-way analysis of variance (one-way ANOVA) followed by LSD comparison test using GraphPad Prism 8 software packages, 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. Additionally, clustering heatmaps were further used to provide an intuitive visualization of selecting candidate metabolites. Finally, all discriminating markers were embedded into a network plot.

Metabolites identification and metabolic pathway analysis

For identification of potential metabolites, the structure information, exact molecular weights and the accurate mass spectrometric fragments with the metabolites were qualified by searching the following freely accessible online biochemical databases: (a) HMDB; (b) METLIN; (c) MassBank; and (d) KEGG.

MetPA (Metabolomics Pathway Analysis), as a comprehensive and user-friendly web-based tool, was used for holistic PATHWAY analysis and mapping of selected metabonomic data.

Results

The overall health status, paw thickness and inflammatory cytokines

Compared with the AS group, the body weight was significantly increased, and the gait score was significantly reduced after 4-week moxibustion (p < 0.05). (Fig. 1A and C). Moxibustion also could significantly reduce the scores of behavioral observation scale compared with AS group (p < 0.01) (Fig. 1D).

Fig. 1

Effects of moxibustion for (A) body weight, (B) hind paw thickness, (C) gait score, (D) behavioral observation scale score, (E) IL-1β, (F) PGE2, (G) IL-6, and (H) TNF-α. AS, ankylosing spondylitis; HC, healthy control; MA, moxibustion at acupuncture points; MNA, moxibustion at non-acupuncture points. (n = 8 mice per group). Data are expressed as means ± standard deviation. ##p < 0.01: HC vs. AS, *p < 0.05: AS vs. MA, **p < 0.01: AS vs. MA, a: p < 0.05: MA vs. MNA.

Compared with the AS group, the paw thickness in the MA group was significantly decreased after 4-week moxibustion therapy (p < 0.01). (Fig. 1B). Moxibustion also significantly decreased the levels of IL-β, PGE2, IL-6 and TNF-α (p < 0.01) in MA group compared with the AS group, while MNA group failed to do so.

Pattern recognition analysis

The representative UHPLC/Q-TOF-MS total ion chromatogram (TIC) profiles are shown in Supplementary C. After PLS-DA processing, clear clustering of the HC, AS, MA and MNA groups in both ESI+ and ESI- modes is shown in Fig. 2, which suggested that the pathological process of AS modeling induced by intraperitoneal injection of PGs and DDA seriously altered normal endogenous metabolic fingerprints of tissue and urine samples in mice. Moreover, variations of tissue and urine metabolic profiling in the MA group was 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.

Fig. 2

PLS-DA score plots of tissue (T) and urine (U) in positive (T1 and U1) and negative ion modes (T2 and U2). AS, ankylosing spondylitis; HC, healthy control; MA, moxibustion at acupuncture points; MNA, moxibustion at non-acupuncture points.

Identification of potential tissue and urine biomarkers and the changing trends among different groups

After pareto scaling, OPLS-DA score plot of tissue datasets showed a remarkable metabolic distinction both when HC vs. AS (Fig. 3A and B) and AS vs. MA (Fig. 3 C and D). The cumulative R2Y and Q2 of the OPLS-DA model in ESI+ and ESI- modes were both above 0.70; suggesting that the models were good to prediction and reliability. Then 27 discriminatory metabolites in both positive and negative ion models were regarded as potential metabolic profiles of PG-induced AS with FDR < 0.05 (HC vs. AS in one-way ANOVA followed by LSD comparison test), FC > 1.33 or <0.77, VIP > 1 and the correlation coefficient |r| > 0.55. Specifically, moxibustion significantly regulated 11 of 27 metabolites in the tissue of mice (FDR <0.05, AS vs. MA in one-way ANOVA followed by LSD comparison test) (Table 1). Of them, the concentrations of taurine, phenyllactic acid, and decanoyl-L-carnitine were up-regulated, while 7-oxocholesterol, adenylsuccinic acid, arachidonoylglycine,1-Stearoyl-sn-glycerol-3-phosphocholine,1-Oleoyl-sn-glycero-3-phosphocholine,1-Myristoyl-sn-glycero-3-phosphocholine, L-threonine and dihydroxyfumarate were down-regulated (Table 1). In addition, 6 metabolites (taurine, adenylsuccinic acid, 1-Stearoyl-sn-glycerol-3-phosphocholine, 1-Oleoyl-sn-glycero-3-phosphocholine,1-Myristoyl-sn-glycero-3-phosphocholine, and phenyllactic acid) were significantly altered in both MA and MNA groups compared with the AS group, most of which were related to the TCA cycle and energy metabolism (Table 1).

Fig. 3

OPLS-DA scores plots for tissue (A, B) and urine (E, F) of AS model group versus HC group in ESI+ and ESI- mode. OPLS-DA scores plots for tissue (C, D), and urine (G, H) of AS model group versus MA group in ESI+ and ESI- mode. The clustering heatmap is used to provide an intuitive visualization of selecting candidate metabolites in tissue (I) and urine (J) samples. AS, ankylosing spondylitis; HC, healthy control; MA, moxibustion at acupuncture points.

Table 1

Potential biomarkers and their metabolic pathways.

Matrix/Ionization mode/No	Identification	Formula	HMDB match	Mass (m/z)	R.T. (min)	VIP	Regulation			Related pathway
							AS*	MA#	MNA#–
Tissue/ESI+
1	Decanoylcarnitine	C17H33NO4	HMDB0000651	368.19	207.10	2.18	↓	↑	–	Lipid metabolism
2	Taurine	C2H7NO3S	HMDB0000251	126.02	279.29	1.80	↓	↑	↑	Taurine and hypotaurine metabolism
3	7-oxocholesterol	C27H44O2	HMDB0000501	401.33	34.27	1.30	↑	↓	–	Lipid metabolism
4	Adenylsuccinic acid	C14H18N5O11P	HMDB0000536	464.31	171.83	1.19	↑	↓	↓	TCA cycle
5	Arachidonoylglycine	C22H35NO3	HMDB0005096	361.27	219.56	1.62	↓	↑	–	Glycine, serine and threonine metabolism
6	1-Stearoyl-sn-glycerol 3-phosphocholine	C10H22NO7PR	METPA0517	568.33	167.43	1.07	↑	↓	↓	Lipid metabolism
7	1-Oleoyl-sn-glycero-3-phosphocholine	C26H52NO7P	HMDB0002815	522.34	172.02	1.20	↑	↓	↓	Lipid metabolism
8	1-Myristoyl-sn-glycero-3-phosphocholine	C9H19NO7PR2	METPA0571	468.30	181.39	2.15	↑	↓	↓	Lipid metabolism
Tissue/ESI-
9	Phenyllactic acid	C9H10O3	HMDB0000779	181.05	175.51	1.79	↑	↓	–	Phenylalanine, tyrosine and tryptophan biosynthesis; Phenylalanine metabolism
10	L-threonine	C4H9NO3	HMDB0000167	118.05	335.19	1.53	↑	↓	–	Glycine, serine and threonine metabolism
11	Dihydroxyfumarate	C4H4O6	HMDB0002050	129.06	46.38	1.83	↑	↓	↓	TCA cycle
Urine/ESI+
12	2-Ethoxyethanol	C4H10O2	HMDB0031213	151.10	64.80	1.04	↑	↓	–	Energy metabolism
13	Valine	C5H11NO2	HMDB0000883	142.09	48.89	2.34	↑	↓	↓	Valine, leucine and isoleucine biosynthesis
14	Glycerol	C3H8O3	HMDB0000131	155.00	257.63	1.32	↑	↓	–	Lipid metabolism
15	Delta-Tocotrienol	C27H40O2	HMDB0030008	460.31	97.07	1.15	↓	↑	–	Tocotorienol biosynthesis
16	Stearoylcarnitine	C25H49NO4	HMDB0000848	504.28	231.16	2.31	↓	↑	–	Lipid metabolism
17	S-citramalic acid	C5H8O5	METPA0309	212.06	241.66	2.23	↑	↓	↓	TCA cycle
18	Isocaproic acid	C6H12O2	HMDB0000689	155.04	70.13	2.20	↓	↑	–	Propanoate metabolism
19	3-Hydroxykynurenine	C10H12N2O4	HMDB0000732	247.07	220.10	1.46	↑	↑	–	Phenylalanine, tyrosine and tryptophan biosynthesis
20	Tryptophan	C11H12N2O2	HMDB0000929	249.07	34.66	1.18	↑	↓	↓	Phenylalanine, tyrosine and tryptophan biosynthesis;Tryptophan metabolism
21	Propanoic acid	C3H6O2	HMDB0000237	170.03	329.54	1.43	↓	↓	–	Propanoate metabolism
22	Ornithine	C5H12N2O2	HMDB0000214	115.09	300.90	1.20	↓	↑	–	Urea cycle
23	Serine	C3H7NO3	HMDB0000187	106.05	398.60	1.41	↑	↓	–	Glycine, serine and threonine metabolism
24	Valeric acid	C5H10O2	HMDB0000892	118.086	366.71	1.69	↑	↓	–	Propanoate metabolism
25	L-phenylalanine	C9H11NO2	HMDB0000159	166.09	241.04	1.82	↑	↓	↓	Phenylalanine metabolism
26	Hippuric acid	C9H9NO3	HMDB0000714	404.19	84.96	1.14	↓	↑	–	Phenylalanine metabolism
Urine/ESI-
27	Purine	C5H4N4	HMDB0001366	152.02	111.65	1.28	↑	↓	–	Purine metabolism
28	cAMP	C10H12N5O6P	HMDB0000058	328.04	267.84	1.91	↑	↓	–	Purine metabolism
29	8R-HETE	C20H32O4	METPA0547	379.25	224.96	2.09	↑	↓	↓	Arachidonic acid metabolism
30	L-alanine	C3H7NO2	HMDB0000161	88.04	394.33	1.27	↑	↓	–	Valine, leucine and isoleucine biosynthesis
31	Acetone	C3H6O	HMDB0001659	111.01	257.50	1.93	↑	↓	–	Lipid metabolism
32	Guanosine	C10H13N5O5	HMDB0000133	282.08	248.85	1.72	↑	↓	–	Purine metabolism
33	Inosine	C10H12N4O5	HMDB0000195	267.07	257.72	1.52	↑	↓	–	Purine metabolism
34	Sebacic acid	C10H18O4	HMDB0000792	201.11	300.81	1.37	↑	↓	–	Lipid metabolism
35	Formylanthranilic acid	C8H7NO3	HMDB0004089	164.04	62.18	1.68	↑	↓	↓	Tryptophan metabolism
36	N-Acetyl-L-glutamate	C7H11NO5	HMDB0001138	188.06	298.48	1.04	↑	↓	↓	Arginine biosynthesis
37	p-Cresol	C7H8O	HMDB0001858	107.05	112.81	1.91	↑	↑	–	Phenylalanine metabolism

AS, untreated ankylosing spondylitis group; HC, healthy control group; MA, moxibustion at acupuncture points; MNA, moxibustion at non-acupuncture points; RT, Retention Time; VIP, Variable importance in the projection.

Differential metabolites: (↓) down-regulated, (↑) up-regulated, and (–) no significant change.

Compared with the HC group, *p < 0.05.

Compared with the AS group, #p < 0.05.

For the urine metabolomics study, OPLS-DA model showed a visible separation between the HC and the AS groups (Fig. 3E and F), and between the AS and MA groups (Fig. 3G and H).

These two OPLS-DA models were both shown to be good fitness and robust (parameters of R2Y and Q2 were both above 0.70). By comparison among different groups, a total of 55 candidate altered variables were screened out as potential metabolites of AS according to the same pre-defined criteria in tissue metabolomics study. In addition, among these 55 differential endogenous metabolites, 26 metabolites were adjusted and contributed to the therapeutic effect of moxibustion, including 7 elevated metabolites (delta-Tocotrienol, isocaproic acid, stearoylcarnitine, 3-hydroxykynurenine, hippuric acid, ornithine and p-Cresol) and 19 decreased metabolites (2-ethoxyethanol, valine, glycerol, S-citramalic acid, tryptophan, propanoic acid, serine, valeric acid, L-phenylalanine, purine, cAMP, 8R-HETE, L-alanine, acetone, guanosine, inosine, formylanthranilic acid, N-Acetyl-L-glutamate, and sebacic acid) in the MA group compared with those in the AS model group (Table 1). While only 6 metabolites valine, S-citramalic acid, tryptophan, L-phenylalanine, 8R-HETE, formylanthranilic acid and N-Acetyl-L-glutamate were reversed by both moxibustion on acupuncture points and non-acupuncture points (Table 1).

Moreover, clustering heat-maps were generated to visually highlight the fluctuations in disturbed tissue (Fig. 3I) and urine (Fig. 3 J) endogenous metabolites and trends of the metabolic concentration difference among different groups. The brown color represents up-regulated metabolic contents and blue color means down-regulated metabolic levels.

Metabolic pathway analysis

In order to identify and visualize the potential target metabolic pathways for moxibustion attenuating AS, 37 endogenous biomarkers listed in Table 1 were imported to MetPA 4.0. It is noteworthy that key influential metabolic pathways “TCA cycle (Raw P = 0.0297)”, “Lipid metabolism (Raw P = 0.0385)”, “Amino Acid metabolism: Aminoacyl-tRNA biosynthesis (Raw P = 4.2443E-4), Alanine, aspartate and glutamate metabolism (Raw P = 0.0039), Nitrogen metabolism (Raw P = 0.0153), Glycine, serine and threonine metabolism (Raw P = 0.0267), Taurine and hypotaurine metabolism (Raw P = 0.0297), Phenylalanine, tyrosine and tryptophan biosynthesis (Raw P = 0.0417), and Valine, leucine and isoleucine biosynthesis (Raw P = 0.0417)”, “Intestinal flora metabolism: Propanoate metabolism (Raw P = 0.0114) and Phenylalanine metabolism (Raw P = 0.0225)”, and “Purine metabolism (Raw P = 0.0350)” significantly contributed to the anti-AS effects of moxibustion (Fig. 4). Finally, histograms of all 37 differential metabolites identified from tissue and urine samples were embedded into the holistic metabolic pathway network graph (Fig. 5).

Fig. 4

Pathway analysis of significant metabolites. a: Aminoacyl-tRNA biosynthesis; b Alanine, aspartate and glutamate metabolism; c: Propanoate metabolism; d: Nitrogen metabolism; e: Phenylalanine metabolism; f: Glycine, serine and threonine metabolism; g: Citrate cycle (TCA cycle); h: Taurine and hypotaurine metabolism; i: Purine metabolism; j: Lipid metabolism; k: Phenylalanine, tyrosine and tryptophan biosynthesis; l: Valine, leucine and isoleucine biosynthesis.

Fig. 5

Summaries of metabolic alterations perturbed by AS modeling and moxibustion with histograms of differential metabolites embedded into the network plot. BCAAs, branched chain amino acids; MCFAs, medium-chain fatty acids; SCFAs, short-chain fatty acid.

Discussion

Our findings clearly indicate that moxibustion can significantly improve the overall health-related phenotypes of AS mice, which include the increasing the body weight and decreasing the gait score and behavioral observation scale score for AS mice. Moreover, moxibustion also decreases the paw thickness and tissue levels of IL-1β, PGE2, IL-6 and TNF-α in AS mice. Furthermore, our results indicate that moxibustion reverse the metabolic dysfunction of AS model mice, mainly involving the TCA cycle, Lipid metabolism, Amino Acid metabolism, Intestinal flora metabolism and Purine metabolism.

The levels of dihydroxyfumarate, adenylsuccinic acid and S-citramalic increase in the PG-induced AS model mice group and it suggests a boosted TCA energy supply to compensate for the shortage of body energy. However, moxibustion decreases levels of dihydroxyfumarate, adenylsuccinic acid and S-citramalic acid and it indicates that moxibustion can contribute to recovering the impairment of the TCA cycle via shifting energy metabolism from glycolysis to aerobic oxidation. Moreover, these metabolic changes were also associated with general improvement in the overall health state as depicted from regaining weight and improving the behavioral observation scale scores.

Acetone and sebaric acid significantly ascended in the PG-induced AS model mice group compared with the HC group, which indicated that a boosted fatty acid β-oxidation was mobilized to make a compensation for the shortage of energy sources caused by TCA energy cycle impairment as we reported earlier. The levels of decanoylcarnitine and stearoylcarnitine in the AS model group are markedly down-regulated compared with the HC group, indicating that fatty acid β-oxidation metabolism is significantly altered and anti-oxidant effect is reduced in the AS model. However, these altered metabolic levels have been normalized following moxibustion intervention, and physical activity improvement and paws with a slight degree of thickness have been observed in the MA group than in the AS model group, which suggest that moxibustion contributes to regulating fatty acid β-oxidation and exerting a defensive capacity against lipid oxidant stress damage in AS mice.

The levels of both cholesterol and phosphocholine metabolites were higher in AS mice than in the controls. The accumulation of cholesterol and phosphocholine metabolites could promote the production of inflammatory cytokines, which lead to the oxidative stress induced cartilage dysfunction in AS. However, moxibustion down-regulates the content levels of cholesterol and phosphocholine compounds markedly, which indicates that moxibustion may ameliorate ectopic fat deposition and protect the cartilage from oxidative damage.

The previous research had shown that Amino Acid metabolism, mainly including tryptophan metabolic pathway, glycine, valine, leucine and isoleucine biosynthesis, and serine and threonine metabolism pathway, was closely linked with the pathogenesis of AS. It was reported that the disturbed tryptophan metabolic pathway in AS was frequently correlated with the immune activation and the inhibition of regulatory T cells. Valine and L-alanine, known as branched chain amino acids (BCAAs), have been implicated to stimulate the production of pro-inflammatory cytokines in rheumatic diseases. The dysregulation of glycine, serine and threonine metabolism was frequently associated with the process of HLA-B27 protein misfolding, the decrease of anti-inflammatory properties and the reducing of immunoglobulin production. However, moxibustion appeared to restore the metabolites and phenotypic indicators correlated with Amino Acid metabolism.

The decreased levels of phenyllactic acid and hippuric acid and increased concentrations of L-phenylalanine and P-Cresol in AS model group imply the abnormal metabolism of intestinal flora and phenylalanine, which played an important role in the pathogenesis of AS. Recently, Ji et al. found that herbal-partitioned moxibustion could protect the integrity of intestinal barrier in rats with ulcerative colitis by reversing microbiota profiles and rebalancing the intestinal microbial environment. In line with previous work conducted by Ji et al, moxibustion can ameliorate the dysbiosis of intestinal flora and the perturbed phenylalanine metabolism pathway in AS mice through regulating metabolites associated with intestinal homeostasis.

Short-chain fatty acid (SCFAs) including isocaproic acid and medium-chain fatty acids (MCFAs) including propanoic acid and valeric acid were preferred microbial metabolites for colonocytes, which played a crucial role in decreasing the intestinal inflammation of AS. Moxibustion can significantly alter the levels of SCFAs and MCFAs, suggesting that moxibustion may exert a defensive capacity against intestinal inflammation and improve the barrier function of the colon in AS mice.

It was reported that purine metabolism was notably altered in adult and pediatric spinal spondyloarthritis, and overproduction of metabolites in purine metabolism could aggravate joint and cartilage tissue injury for AS patients. Previously, heat-reinforcing acupuncture had exerted a positive effect on decreasing the purine metabolism on arthritis rabbit with cold syndrome. Results of this study also exhibited that the urine levels of purine, cAMP, guanosine and inosine in AS mice were markedly down-regulated by moxibustion therapy, suggesting a similarity of anti-AS mechanisms to acupuncture.

There are several potential limits of this study. Firstly, moxibustion can ameliorate the dysbiosis of AS mice via regulating the urine metabolites involved in the intestinal flora metabolism. However, metabolites from fecal extracts and intestinal tissues may directly reveal the crosstalk between gut microbiota and intestinal cells. Secondly, several important lipid molecules, including sphingomyelins, neutral lipids, and triglycerides, yet were not detected in this research due to the limitation of LC/MS technique.

In conclusion, moxibustion obviously exerted a reverse effect on AS-induced metabolic alterations, especially the expression of metabolic components involved in TCA cycle, Lipid metabolism, Amino Acid metabolism, Intestinal flora metabolism and Purine metabolism. Thus, the UHPLC-Q-TOF/MS based metabolomics approach, as a novel and powerful tool, can help us to gain the insights into potential mechanisms of action of moxibustion for AS.

Acknowledgement

We would like to appreciate the editor and all staff working in editorial office of IMR. We also appreciate all anonymous reviewers who provided insightful suggestions for our manuscript.

Author contributions

Xiao Xu: Conceptualization, Data curation, Formal analysis, Funding acquisition, Methodology, Resources, Software, Validation, Visualization, Writing - original draft, and Writing - review & editing. Ya-Nan Shi: Conceptualization, Data curation, Formal analysis, Software, Validation, Visualization, Writing - original draft. Rong-Yun Wang: Data curation, Formal analysis, Software, Validation, Visualization. Cai-Yan Ding: Resources, Software, Validation, Visualization and Visualization. Xiao Zhou: Methodology, Resources. Yu-Fei Zhang: Methodology, Resources. Zhi-Ling Sun: Validation, Visualization, Funding acquisition. Zhi-Qin Sun: Conceptualization, Investigation, Methodology, Project administration, Resources, Supervision, Writing - review & editing. Qiu-Hua Sun: Conceptualization, Investigation, Methodology, Project administration, Resources, Supervision, Funding acquisition, Writing - review & editing.

Conflict of interest

The authors declare that there is no conflict of interest.

Funding

This study was supported by National Natural Science Foundation of China (No. 81904274, 81973756 and 81774383).

Ethical statement

This study has been approved by the ethical review of experimental animal welfare in Zhejiang Chinese Medical University (No.SCXK2016-0010).

Data availability

The data will be made available upon reasonable request.