Correlation Analysis between Gut Microbiota and Metabolites in Children with Systemic Lupus Erythematosus

Systemic lupus erythematosus (SLE) is an autoimmune-mediated diffuse connective tissue disease characterized by immune inflammation with an unclear aetiology and pathogenesis. This work profiled the intestinal flora and faecal metabolome of patients with SLE using 16S RNA sequencing and gas chromatography-mass spectrometry (GC-MS). We identified unchanged alpha diversity and partially altered beta diversity of the intestinal flora. Another important finding was the increase in Proteobacteria and Enterobacteriales and the decrease in Ruminococcaceae among SLE patients. For metabolites, amino acids and short-chain fatty acids were enriched when long-chain fatty acids were downregulated in SLE faecal samples. KEGG analysis showed the significance of the protein digestion and absorption pathway, and association analysis revealed the key role of 3-phenylpropanoic acid and Sphingomonas. Sphingomonas were reported to be less abundant in healthy periodontal sites of SLE patients than in those of HCs, indicating transmission of oral species to the gut. This study contributes to the understanding of the pathogenesis of SLE disease from the perspective of intestinal microorganisms, explains the pathogenesis of SLE, and serves as a basis for exploring potential treatments for the disease.


Introduction
Systemic lupus erythematosus (SLE) is an autoimmunemediated diffuse connective tissue disease characterized by immune inflammation [1]. The presence of antinuclear antibodies in the serum and multisystem involvement are the two main clinical features of SLE, and the incidence rate of SLE in women is significantly higher than that in men [2]. Lupus nephritis (LN) occurs when SLE is complicated by renal damage. Approximately 60% of SLE patients suffer from LN complications that cause a higher mortality than SLE patients without LN [3]. At present, the aetiology and pathogenesis of SLE and LN are still unclear [4]. Studies have shown that the occurrence and development of some autoimmune diseases, such as inflammatory bowel disease [5,6], rheumatoid arthritis [7,8], and SLE [9,10], are correlated with changes in intestinal microecology. Some studies found that the intestinal microbes of patients or animals with SLE were different from those of the control group through experimental models and clinical models [11][12][13], and the reconstruction of intestinal microecology through diet and drug regulation could help improve the severity of the disease [14]. There are many hypotheses on the interaction mechanism between intestinal microbes and the host and its relationship with SLE, but there are no mature theories yet. Currently, an increasing number of studies have shown that metabolites of the intestinal flora, such as short-chain fatty acids (SCFAs), free fatty acids (FFAs), amino acids, and arachidonic acid, correlate with autoimmune reactions [14][15][16]. However, the relationship between microbial metabolites and SLE remains unknown, and previously published studies are limited to adult patients. In this study, we collected the intestinal flora and metabolite data of 33 SLE children and 28 healthy controls (HCs) to clarify the pathogenesis of the disease from the perspective of intestinal microorganisms as much as possible, explain the pathogenesis of SLE, and serve as a basis for exploring potential treatments for the disease.

Material and Methods
2.1. Ethical Statement. All enrolled patients and healthy controls understood the sampling process and study plan in detail before sampling and signed informed consent forms. Ethical approval for this study was obtained from the Ethics Committee of the Second Xiangya Hospital of Central South University.

Study Subjects.
A total of 28 healthy controls and 33 SLE patients were enrolled from the Second Xiangya Hospital of Central South University between December 2018 and April 2019. The diagnostic criteria of SLE were consistent with SLE classification criteria 2 revised by American College Rheumatology (ACR) in 1997 [17]. SLE patients were given hormone and/or immunosuppressive therapy. None of the selected subjects had applied or taken antibiotics, probiotics, symbiotic, or other microecological preparations for at least 3 months before sampling.
2.3. 16S rRNA Microbiota Analysis. Total genomic DNA from samples was extracted using the CTAB/SDS method. 16S/18S rRNA genes were amplified using primers including 16S V4-V5 : 515F-907R, 18S V9 : 1380F-1510R, and ITS1 : ITS1F-ITS2R. All PCRs were carried out in 30 μl reactions with 15 μl of Phusion® High-Fidelity PCR Master Mix (New England Biolabs), 0.2 μM forward and reverse primers, and approximately 10 ng of template DNA. Thermal cycling consisted of initial denaturation at 98°C for 1 min, 30 cycles of denaturation at 98°C for 10 s, annealing at 50°C for 30 s, and elongation at 72°C for 60 s then 5 min. Electrophoresis was performed on a 2% agarose gel. PCR products were purified with GeneJET Gel Extraction Kit (Thermo Scientific). Sequencing libraries were generated using the NEBNext® Ultra™ DNA Library Prep Kit for Illumina (NEB, USA), and the library quality was assessed on the Qubit@ 2.0 Fluorometer (Thermo Scientific) and Agilent Bioanalyzer 2100 system. The library was sequenced on an Illumina MiSeq platform.
A total of 5507169 sequences were used for analysis. Sequence analysis was performed by the UPARSE software package using the UPARSE-OTU and UPARSE-OTUref algorithms. In-house Perl scripts were used to analyse the alpha diversity and beta diversity. Cluster analysis was preceded by principal component analysis (PCA). QIIME calculates both weighted and unweighted UniFrac distances. We used unweighted UniFrac distance for principal coordinate analysis (PCoA) and unweighted pair group method with arithmetic mean (UPGMA) clustering. Metastats software was utilized to confirm differences in the abundances of individual taxa between the two groups, and LEfSe was used for the quantitative analysis of biomarkers within different groups. To identify differences in microbial communities between the two groups, ANOSIM [18] and MRPP (multiresponse permutation procedure) [19] were performed based on Bray-Curtis dissimilarity distance matrices.

GC-MS Analysis.
The instrumental status and balance as well as the stability of the gas chromatography-mass spectrometry system were evaluated using QC samples. Samples (60 mg) were mixed with 200 μg water and vortexed for 60 s before being added to 800 μl methanol acetonitrile (Merck, 1499230-935) solution (1 : 1, v/v), followed by vortexing for 60 s and low-temperature ultrasound for 30 min twice. Samples were stored at -20°C for 1 h to precipitate protein.
Finally, the samples were centrifuged for 20 min at 14000 RCF at 4°C, and the supernatant liquor was collected for freeze drying and stored at -80°C.
For GC analysis, the samples were separated by an Agilent 1290 Infinity LC ultra-high-performance liquid chromatography (UHPLC) HILIC column at 25°C with a flow rate of 0.3 ml/min. The mobile phase was the composite of A (water, 25 mm ammonium acetate, 25 mm ammonium) and B (acetonitrile). Gradient elution procedures were as follows: B was 95% in 0-1 min, then changed linearly from 95% to 65% in 1-14 min and 65% to 40% in 14-16 min, and maintained at 40% in 16-18 min before changing linearly from 40% to 95% in 18-18.1 min, maintaining at 95% in 18.1-23 min. The samples were placed in an automatic sampler at 4°C throughout the analysis. A random sequence was used for the continuous analysis. QC samples were inserted into the sample queue for monitoring and evaluation of the system.
For Q-TOF MS analysis, electrospray ionized positive ions and negative ions were detected. The samples were separated by UHPLC and analysed by a Triple TOF 5600 mass spectrometer (AB SCIEX). The ESI source conditions after HILIC chromatographic separation were as follows: ion source gas 1 (Gas1): 60, ion source gas 2 (Gas2): 60, curtain gas (CUR): 30, source temperature: 600°C, ion spray voltage floating (ISVF) ±5500 V (positive and negative modes), TOF MS scan m/zrange: 60-1000 Da, product ion scan m/z range: 25-1000 Da, TOF MS scan accumulation time 0.20 s/spectra, and product ion scan accumulation time 0.05 s/spectra. The second-order mass spectrum was obtained by information-dependent acquisition (IDA) with a high sensitivity mode, a declustering potential (DP): of ± 60 V (positive and negative modes), collision energy of 35 ± 15 eV, and IDA excluded isotopes within 4 Da with 6 candidate ions to monitor per cycle.
The original data were converted into the mzXML format by ProteoWizard, and then, the XCMS program was used for peak alignment, retention time correction, and peak area extraction. Accurate mass number matching (<25 ppm) and secondary spectrogram matching were used to identify the metabolite structure, and the laboratory database was retrieved. After pretreatment by Pareto scaling, multidimensional statistical analysis was performed, including principal component analysis (PCA), supervised partial least square discriminant analysis (PLS-DA), and orthogonal partial least square discriminant analysis (OPLS-DA). One-dimensional statistical analysis included Student's t-test and multiple analysis of variation.

Association
Analysis. The correlation analysis is mainly based on statistical algorithms to find the correlation between significantly different flora (p < 0:05) and the significantly different metabolites obtained through nontarget metabolomics analysis. Principal component analysis (PCA) and multivariate statistical analysis were performed with SIMCA version 14.1. KEGG pathway analysis was based on the online Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.kegg.jp/). Differentially abundant proteins/modified peptides/genes/and metabolites/lipids were log 2 scaled (TMT/iTRAQ) or Z-score scaled (label free) and concatenated into one matrix. Then, the correlation coefficients among all the molecules in the matrix were calculated with the Pearson algorithm in R version 3.5.1. Cytoscape version 3.5.1 was used to calculate the correlation network.  Table 1, including their sex, age, disease activity parameters, autoantibody status, and treatment. According to two-tailed unpaired Student's t-test, BMI was not significantly different between the SLE group and HC group (p = 0:234). All the children in the SLE group were also diagnosed with LN.

Intestinal Flora Is Different between SLE Patients and
Healthy Groups. After data optimization, the average length of 99.87% of the sequence bars was within the range of 401-450 bp, and a total of 2,159 operational taxonomic units (OTUs) were retained for subsequent analysis.
Using a t-test, we compared the Chao1 index (p = 0:327) of the two groups and found no significant difference, suggesting that the species abundance does not change in SLE patients. Additionally, Shannon (p = 0:214), coverage (p = 0:988), and Simpson (p = 0:414) indexes that estimate community diversity also showed similar results. Taken  Figure 1: Bray-Curtis statistical algorithm was used to calculate the distance between each sample and obtain the distance matrix. The distance between the SLE group and the HC group with quartiles was calculated, and a box diagram of the distance between the two groups was drawn to compare the distance distribution differences. Student's two-sample t-test method was used to compare the distance between each group.
3 Journal of Immunology Research together, these results show that the richness of the intestinal flora did not change among SLE patients.
As we can see in Figure 1, the sample distance of the control vs. control was significantly shorter than SLE vs. SLE group and control vs. SLE group (p < 0:05). Interestingly, although the alpha diversity of SLE patients remained unchanged compared with that of healthy controls, the sample distance of the former was obviously longer, suggesting healthy children possess microbiota with more similar species composition and the heterogeneity of the microbiota among SLE patients.
A scatter plot based on principal coordinates analysis (PCoA) revealed a clear separation between the SLE and HC groups ( Figure 2). It should be highlighted that a principal component analysis (PCA) of SLE and HC subjects, based on 16 rRNA microbiota profiles, did not show distinct clustering patterns. We witnessed the association between gastrointestinal microbiota and immune factors at the level of bacterial composition but not at the level of the metabolite landscape of gut microbiota, which suggests that SLE can influence the heterogeneous species inhabiting the gut but not the metabolic pattern [20].
Using LEfSe software, a nonparametric factorial Kruskal-Wallis (KW) sum-rank test was conducted to identify species with distinct abundance differences before undertaking linear discriminant analysis (LDA) to estimate the influence of each component (species) on the difference. The difference in microflora from phylum to genus was determined by an LDA value greater than 2 and a p value less than 0.05 (Table 2). We can see that for SLE patients, Proteobacteria at the phylum level, Gammaproteobacteria at the class level, Enterobacteriales at the order level, and Enterobacteriaceae at the family level are the most relevant components, and they possess a subordinate relationship. On the other hand, Ruminococcaceae, Christensenellaceae, and Family_XIII at    Journal of Immunology Research  (Figure 3).
The nonparametric factorial Kruskal-Wallis sum-rank test was used for calculating the difference of abundance. Log value of abundance reflects the abundance of the different species. LDA value estimates the influence of the abundance of each component (species) on the overall abundance differences. p < 0:05 and LDA > 2 represent significantly different species.

Different Metabolites' Analysis by GC-MS.
Variable importance for the projection (VIP) obtained from the OPLS-DA model was used to measure the influence intensity and explanatory ability of the expression pattern of each metabolite on the classification and discrimination of samples of each group, and the different metabolites with biological significance were excavated.
This study defined metabolites with VIP > 1 in the multidimensional statistical analysis and p < 0:05 in univariate statistical analysis as significantly different metabolites, including 24 in positive ion mode and 37 in anion mode. Ruling out 6 overlapping data, there were a total of 55 significantly different metabolites (Table 3), including mainly long-chain fatty acids and amino acids. Among them, cyclohexylsulfamate had the highest VIP value, indicating that its decline had a greater impact on the occurrence of SLE diseases.
Of the 55 SLE-altered metabolites submitted to the KEGG (Kyoto Encyclopedia of Genes and Genomes) website (http://www.kegg.jp/), we were able to match 27 of them with 66 pathways. The p value represents the significance of enrichment of the KEGG pathway, and the number of differentially expressed metabolites contained in the KEGG pathway reflects the degree of influence of SLE on each pathway to some extent (Figure 4). Taking two factors into consideration, protein digestion and absorption (p = 1:75E − 9) were the most significant among all these pathways, involving Ltryptophan, tyramine, L-phenylalanine, L-leucine, L-methionine, L-alanine, L-glutamine, L-valine, L-isoleucine, and Ltyrosine.  Figure 3: Differential analysis using LEfSe software. Red nodes represent the microbial group that plays an important role in the HC group, green nodes represent the microbial group that plays an important role in the SLE group, and yellow nodes represent the microbial group that does not play an important role in both groups. 6 Journal of Immunology Research The network of metabolites and intestinal flora using Cytoscape version 3.5.1 is shown in Figure 5. We can see that the most central location belongs to 3-phenylpropanoic acid, which exhibits a strong positive correlation with Family XIII UCG-001 and a negative correlation with all 5 differentially expressed genera of Proteobacteria. Sphingomonas of the Proteobacteria phylum plays the most important role as a bacterium, and it is obviously positively associated with all 10 metabolites in the protein digestion and absorption pathway mentioned above.

Discussion
The current study found that alpha diversity of the SLE group is not different when compared to healthy controls, and the beta diversity was partially altered, which was in accordance with a previous study among adult SLE patients [20]. The fact that the alpha diversity is not statistically different between the healthy control group and the group affected by the pathology could be related to age-dependent factors and also to the small number of children in the analysis. To be precise, PCoA successfully identified 2 separated clusters of SLE and HC subjects when PCA did not, indicating that the changes induced by SLE became marked at the functional level, i.e., the level of microbial population structure (or 16S rRNA), regardless of the heterogeneities that appeared above at the highest level of the functional hierarchy, i.e., the metabolite level ( Figure 1). In contrast to our outcomes on the alpha and beta diversity, the outcomes of a study by Zhu revealed upregulated alpha diversity in the SLE group; in addition, Rojo et al. reported beta diversity in PCA but not PCoA [16,21]. This inconsistency may be due to age and sex differences since previous papers contained only adult patients and fewer or no male patients.
Another important finding was the increase in Proteobacteria and Enterobacteriales at the phylum and family levels among SLE patients in our project in southern China, and this finding was also reported by Wei et al., He et al., and Hevia and Milani from northern China. However, results from studies conducted in Spain and southern China do not support our results [22][23][24]. In Spain, the intestinal flora of SLE patients was characterized by a significantly lower Firmicutes/Bacteroidetes ratio and showed a depletion of Lachnospiraceae and Ruminococcaceae and an enrichment of

Journal of Immunology Research
Bacteroidaceae and Prevotellaceae [24]. In southern China, the intestinal flora of SLE patients was characterized by a significant increase in Actinobacteria [22]. It can therefore be assumed that the flora composition is irrelevant to the geographical location. On the other hand, consistent with the literature, this research found that Ruminococcaceae levels were higher in healthy controls [21,23,24]. Moreover, there was a list of other genus assigned to the Lachnospiraceae family that are widely descripted as beneficial bacteria in literature elevated in health controls, including Lachnospira-ceae_UCG_008 and Lachnospiraceae_UCG_004 [25,26]. Additionally, previous research has established that the increase in Proteobacteria and the decrease in Ruminococcaceae were associated with lupus nephritis with gastrointestinal damage [27]. Na et al. proposed that an increased prevalence of the bacterial phylum Proteobacteria is a marker for an unstable microbial community (dysbiosis) and a potential diagnostic criterion for disease [28]. Shang et al. identified Ruminococcaceae to be responsible for the degradation of diverse polysaccharides and fibres [29].
An obvious finding to emerge from the analysis is that amino acids were enriched in SLE faecal samples compared to HCs (Table 3), which was in accordance with previous research that found that their plasma concentration decreased [30][31][32][33]. Using gas chromatography-mass spectrometry, Yan et al. were able to detect the serum depletion of 12 amino acids, including tryptophan, alanine, proline, glycine, serine, threonine, aspartate, glutamine, asparagine, lysine, histidine, and tyrosine [30]. Similarly, Ouyang et al. found that L-alanine, L-isoleucine, L-lysine, L-phenylalanine, L-tyrosine, and L-valine were downregulated in the plasma of SLE patients [31]. Notably, four essential amino acids, Lvaline, L-leucine, L-phenylalanine, and L-tryptophan, were listed among the metabolites depleted in plasma but enriched in faeces. These four amino acids have important physiological functions, such as the regulation of immunity, metabolism, and neural activity. Choi et al. reported that low dietary tryptophan prevented autoimmune pathology in lupus-prone mice, whereas high dietary tryptophan exacerbated disease [32]. Moreover, glucogenic amino acids, such as proline and L-methionine, and glucogenic and ketogenic amino acids, such as L-tyrosine, were increased in the faeces of SLE patients, indicating that there might be disorders in glucose metabolism and energy metabolism, as these amino acids could emerge as potential energy sources [33,34]. Regarding the fatty acids, we found the abundance varied, as most long-chain fatty acids decreased while short-chain fatty acids increased in SLE children. This finding is contrary to a recent study that reported enriched fatty acids in the faeces of the SLE group [35]. Another study also established that long-chain fatty acids were less abundant in SLE patient serum [36]. At least half of SLE patients have gastrointestinal  Figure 5: Circles represent significantly different genera, and rectangles represent significantly different metabolites. The color of the line represents the positive and negative value of correlation coefficient (blue represents a negative correlation and red represents a positive correlation), and the thickness of the line is proportional to the absolute value of the correlation coefficient. The node size is positively correlated with its degree (the greater the degree, the larger the node size). 9 Journal of Immunology Research symptoms, including nausea, vomiting, anorexia, abdominal pain, diarrhoea, and abdominal distension, which may lead to an imbalance of amino acids and fatty acids [37].
The results of this study show that SLE-enriched faecal amino acids were significantly located in the protein digestion and absorption, ABC transporters, and aminoacyl-tRNA biosynthesis pathways. Amino acids are mainly absorbed in human intestinal mucosal cells through carrier proteins and the γ-glutamine cycle. Through type I ABC transporters in the third subgroup, intestinal prokaryotes and archaea obtain amino acids, which can be used for the synthesis of microbial proteins, including the biosynthesis of aminoacyl-tRNA, microbial energy metabolism, and conversion into a variety of physiologically active substances, such as neurotransmitters and hormones [38]. As a result, the increase in amino acids in faeces can impact the physiological activities of microorganisms and their hosts. For example, SLE-rich faecal branched-chain amino acids, including L-leucine and L-valine, can be used as precursors of membrane fatty acids, as well as key synergistic regulators of pathogen growth and virulence. Low-proline or lowprotein diets in germ-free mice parasitized by human intestinal flora led to reduced expansion of wild-type C. difficile after challenge, indicating that the effectiveness of amino acids may be important for C. difficile infection [39]. Tryptophan is essential for the growth of certain pathogens, and severe tryptophan deficiency may hinder normal Chlamydia trachomatis onset and reduce its activation [40].
Of association analysis, this study found that Sphingomonas plays the most important role in the network of metabolites and intestinal flora. We found higher levels of the bacterium in SLE children. It is encouraging to compare this rise with that found by Corrêa et al., who found that Sphingomonas was at lower relative levels in healthy periodontal sites of SLE patients [41]. Chen et al.'s research provided evidence for increased microbial transmission along the gastrointestinal tract of most certain species in SLE patients compared to HCs, implicating a damaged barrier in the upper gastrointestinal tract in SLE patients who possibly cause transmissible oral species to colonize the gut [42].
It is unfortunate that most studies concerning gut microbiota and metabolites in SLE and LN patients include a study cohort of only adults, so our literature search may not be comprehensive enough. There is also a potential confounding effect of medications, as most patients were receiving treatment at the time of sample collection, and a large list of medications is postulated to affect the balance of taxa within the gut microbiome community [43,44]. Another limitation of this study is that the SLE group and the HC group were not perfectly matched.

Conclusion
This study set out to investigate intestinal flora and metabolomics in SLE children. This study has shown that Proteobacteria was increased in SLE patients. Of this phylum, the genus Sphingomonas plays a very important role, as it exerts a significant positive correlation with a series of proteins in the protein digestion and absorption pathway. Sphingomonas were reported to be less abundant in healthy periodontal sites of SLE patients than in those of HCs, indicating the transmission of oral species to the gut. Adjunct therapeutic modalities of adjusting microbiota can be innocuous and worth trying. Considerably, more work will need to be done to determine the pathogenesis and treatment of SLE in children.

Data Availability
All data are available in our manuscript.

Ethical Approval
Ethical approval for this study was obtained from Ethics Committee of the Second Xiangya Hospital of Central South University.