16S rRNA Gene Sequencing of Gut Microbiota in Rheumatoid Arthritis Treated with 99Tc-MDP

Introduction

Rheumatoid arthritis (RA) is a systemic autoimmune disease. As the disease progresses, it will cause the gradual loss of joint function, which will seriously affect the life and work of patients. The global prevalence of RA is about 1%, with a much higher incidence in women than in men.¹ Although the pathogenesis of RA is unclear, scholars agree that it is closely related to environmental factors, genetic susceptibility, and immune dysfunction. Gut microbes happen to interact closely with all the three.^2–4 Lifestyle, dietary preferences, and gut flora metabolites affect the development and prognosis of RA.^3,5 Meanwhile, the RA susceptibility gene HLA-DRB1 favors a more inflammatory gut microbiome.⁶ In addition, gut microbes are indispensable for the maturation of the immune system, especially mucosa-associated lymphoid tissue.⁷ In recent years, with the development of high-throughput sequencing technologies, more and more evidence has shown that the pathogenesis of RA is closely related to the host gut microbiota, especially in the induction of chronic inflammation and activation of autoimmunity.^8,9

Technetium methylene diphosphonate (99Tc-MDP) is a nuclide preparation, and its main component is the complex formed by the trace element technetium and methylene diphosphonic acid (Indications: RA and orthopaedic diseases). 99Tc-MDP has strong anti-inflammatory and bone repair effects, and at the same time has high safety. It is often used in the treatment of RA in China.^10–12 99Tc-MDP can also reduce the secretion of inflammatory cytokines such as TNF-α and IL-6 by regulating immune cells.^13,14 Clinically, 99Tc-MDP is often used in combination with methotrexate (MTX), the basic drug for RA. The results of a multicenter, double-blind RCT showed that 99Tc-MDP and MTX combination therapy was more effective in RA patients and did not increase the incidence of adverse events compared with 99Tc-MDP or MTX monotherapy.¹⁵ There have been some studies on MTX and gut microbial interactions, but there are no studies on the effect of 99Tc-MDP on the gut flora of RA patients.^16,17 Therefore, we investigated the effect of 99Tc-MDP on the gut microbiota using 16S rRNA gene sequencing and performed bioinformatics analysis. This study provides new information on the mechanism of 99Tc-MDP in RA treatment.

Materials and Methods

Participants

All inpatients who were definitively diagnosed with RA based on clinical, laboratory and imaging status from September to November 2019 at the Department of Rheumatology and Immunology, the Second Affiliated Hospital of Guizhou University of Traditional Chinese Medicine were screened. The RA diagnoses also fulfilled the 2010 American College of Rheumatology/ European League Against Rheumatism (ACR/EULAR) classification criteria for RA.¹⁸ The age, tender joint count (TJC), swollen joint count (SJC), erythrocyte sedimentation rate (ESR), C-reactive protein (CRP), rheumatoid factor (RF), anti-cyclic citrullinated peptide antibody (ACPA), and other clinical conditions of all the participants were recorded. The patients were advised to maintain their usual diet, lifestyle, and exercise throughout the study period. The study was approved by the Ethics Committee of the Second Affiliated Hospital of Guizhou University of Traditional Chinese Medicine (GK2019004), and all participants signed an informed consent form. All study procedures involving human participants were in accordance with the 1964 Declaration of Helsinki and its subsequent revisions.

The patients were enrolled according to the following criteria: (1) Lived in Guizhou Province in the previous year. (2) Aged less than 75 years old and more than 18 years old, and were able to take care of themselves. (3) Oral MTX (10 mg/week) monotherapy for more than 2 months. (4) Disease activity score DAS28-ESR > 3.2.¹⁹ (5) After evaluation by two senior rheumatologists, 99Tc-MDP could be added to the treatment and there was no drug contraindication. The patients were excluded based on any of the following factors: (1) History of severe gastrointestinal diseases such as gastroenteritis, peptic ulcer, and gastric tumor. (2) Antibiotics, glucocorticoids, disease-modifying antirheumatic drugs (DMARDs) except MTX and other immunosuppressive drugs were used within 2 months. (3) Special diet within 2 weeks, such as Chinese herbal medicine, probiotics, yogurt, coffee, and wine. (4) Complicated with hypertension, hyperlipidemia, diabetes, and other serious cardiovascular diseases or metabolic diseases. (5) Combined with other autoimmune diseases such as systemic lupus erythematosus, ankylosing spondylitis, and multiple sclerosis. (6) Children, pregnant or lactating women, and patients with mental illness. (7) Those who discontinued medication without authorization or were lost to follow-up.

Sample Collection

Fresh faecal samples were collected from RA patients on the first day of admission and transported by cold chain to a freezer at −80°C for storage within 10 minutes. These samples were labelled as before the 99Tc-MDP treatment group (Before). After 10 days of treatment with the combination of 99Tc-MDP and MTX, fresh stool samples were collected again using the same method and labelled as after the 99Tc-MDP treatment group (After).

Details of treatment for RA patients: (1) 99Tc-MDP: agent A (5 mL, containing technetium [⁹⁹Tc] 0.05 μg) and agent B (containing methylene diphosphonic acid of 5 mg and stannous chloride of 0.5 mg, 3 vials, a total of 16.5 mg) fully mixed and then set for 5 minutes. Added to 250 mL of 0.9% sodium chloride solution, intravenous drip, once a day for 10 consecutive days. (2) MTX: Patients had been treated with MTX for more than 2 months before hospitalization; therefore the dosage remained unchanged according to the previous dosage (10 mg/week orally).

DNA Extraction, Amplification, Library Construction and Sequencing

Total genomic DNA from samples was extracted using the TIANamp Stool DNA kit (TLANGEN Co., China). DNA concentration and integrity were measured with a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, USA) and agarose gel electrophoresis, respectively. PCR amplification of the V3-V4 hypervariable regions of the bacterial 16S rRNA gene using universal primer pairs (343F: 5′-TACGGRAGGCAGCAG-3′; 798R:5′-AGGGTATCTAATCCT-3′). The reverse primer contained a sample barcode and both primers were connected with an Illumina sequencing adapter.

The PCR products were purified with Agencourt AMPure XP beads (Beckman Coulter Co., USA) and quantified using Qubit dsDNA assay kit. Sequencing was performed on an IlluminaNovaSeq6000 with two paired-end read cycles of 250 bases each. (Illumina Inc., San Diego, CA, USA; OEBiotech Company, Shanghai, China). Paired-end reads were filtered low-quality sequences, denoised, merged, and detect and cut off the chimera reads using DADA2 with the default parameters of QIIME2.^20,21 Last, the software output the representative reads and the ASV abundance table. All representative reads were annotated and blasted against the Silva database (Version 138) using q2-feature-classifier with the default parameters.

Bioinformatics Analysis

Differences between the two groups were also counted at the levels of phylum, class, order, family, and genus. To evaluate the 99Tc-MDP treatment infection on the diversity of gut microbiota, the alpha-diversity that includes Chao1 index and Shannon index was calculated.²² The similarity of bacterial communities was estimated using beta diversity. The Bray-Curtis dissimilarity and Weighted UniFrac dissimilarity were calculated.

LEfSe analysis, linear discriminant analysis (LDA) coupled with effect size measurements, was used to identify the main differential bacteria. LDA score indicates the magnitude of the difference contribution, the larger the absolute value the greater the difference contribution. Random Forest is a machine learning classification model that can efficiently and accurately classify microbial community samples. Mean Decrease Gini is used to compare the importance of variables by calculating their effect on the heterogeneity of the classification model. A higher value indicates that the variable (species) is more important.

Association analysis was performed on the top 30 genera with relative abundance in the gut microbiota to determine the correlation between species (Spearman correlation coefficient). KEGG pathway enrichment analysis was performed at three levels. The Clusters of Orthologous Groups of proteins (COG) database was used to infer the main biological functions of the microbiota. The brief process is shown in Figure 1.

Figure 1 Brief flowchart.

Statistical Analysis

Statistical analyses of clinical status of the RA patients were performed using IBM SPSS Statistics for Windows, version 25 (IBM Corp., Armonk, N.Y., USA). Continuous variables with a normal distribution and homogeneity of variance are represented by mean ± standard deviation (SD), and comparisons between groups were performed by t-test. Continuous variables that do not conform to the normal distribution are represented by the median (interquartile range) and comparisons between groups were performed by the Mann–Whitney U-test. Categorical variables are expressed in terms of quantity and percentage (n, %). The above test methods were all statistically significant with P < 0.05.

Results

Clinical Status

A total of 64 fresh stool samples were collected from 32 RA patients. After the combined treatment of 99Tc-MDP and MTX, the ESR, CPR, RF, and DAS28-ESR of RA patients were significantly decreased (P < 0.05). It indicated that the patients had reduced disease activity and improved disease severity after treatment (Table 1).

Table 1 Clinical Status Before and After Treatment

Sequencing Data

We sequenced a total of 64 stool samples and obtained 3,198,961 clean tags after quality control. After removing chimeras, 2,888,964 valid tags (90.3%) were obtained. The valid tags of each sample were from 29,578 to 49,022. The mean length of the valid tags was from 409 bp to 424 bp. The number of OTUs of each sample was from 331 to 1680. Annotated information based on matching OTUs included 19 phyla, 47 classes, 101 orders, 175 families, 413 genera, and 138 species of gut microbes.

99Tc-MDP Therapy Regulated the Abundance of Gut Microbiota

The gut microbial composition of RA patients before and after 99Tc-MDP treatment is shown in Figure 2. At the phylum level (Figure 2a), Bacteroidetes, Firmicutes, Proteobacteria, Actinobacteria, and Fusobacteria were the five most abundant phyla, accounting for over 98% of all the microorganisms. Of these, Bacteroidetes abundance increased while Firmicutes abundance decreased. At the genus level (Figure 2b), Bacteroides, Escherichia-Shigella, Prevotella_9, Faecalibacterium, and Streptococcus were the five most abundant genera. Of these, Escherichia-Shigella decreased by 30% and Prevotella_9 increased by 55%.

Figure 2 Compositional of gut microbes. (a) Relative abundance at the phylum level. (b) Relative abundance at the genus level.

Compared to the Before group, the Deferriactors increased significantly at the phylum level (Figure 3a). At the class level, the abundance of Deferribacteres increased significantly, while Erysipelotrichia decreased (Figure 3b). At the order level, Deferribacterales and Solibacterales increased, while Erysipelotrichales, Actinomycetales, and Acetobacterales decreased (Figure 3c). At the family level, Enterococcaceae, Deferribacterales, Solibacteraceae_(Subgroup_3), Promicromonosporaceae, and Bacteroidales_RF16_group increased, while Erysipelotrichales, Saprospiraceae, and Acetobacterales decreased (Figure 3d). At the genus level, Acinetobacter, Enterococcus, and Prevotella_6 increased, whereas Ruminococcaceae_UCG-014, Lachnospiraceae_UCG-010, Butyricicoccus, Coprococcus_2, Hungatella, Actinomyces, and Anaerovibrio decreased (Figure 3e).

Figure 3 Statistically different species (Top 10). (a) Differential species at the phylum level. (b) Differential species at the class level. (c) Differential species at the order level. (d) Differential species at the family level. (e) Differential species at the genus level.

Gut Microbial Diversity Unaffected with 99Tc-MDP Therapy

The alpha-diversity results showed that the Chao index and Shannon index had no significant differences between the two groups (P > 0.05; Figures 4a and b). The similarity of bacterial communities was estimated using beta diversity. Similarity analysis showed that there was no significant difference between the two groups (Bray-Curtis, R2=0.0129, P > 0.05; Weighted UniFrac, R2=0.0107, P > 0.05). Principal Coordinates Analysis (PCoA) was performed based on Bray–Curtis distances and weighted UniFrac distances to visualize the similarity of the bacterial community structure in the two groups. As shown in Figures 4c and d, the clustering and bacterial composition of the two groups were similar.

Figure 4 Biodiversity analysis. (a) Chao1 index of α-diversity. (b) Shannon index of α-diversity. (c) Bray-Curtis dissimilarity (R2=0.0129, P > 0.05). (d) Weighted UniFrac dissimilarity (R2=0.0107, P > 0.05).

LEfSe Analysis to Identify Main Differential Bacteria

LEfSe analysis was used to reveal the composition of species that differed between two groups of biomes. It included the magnitude of the species’ contribution to the difference (Figure 5a) and cladogram annotation of the differential species (Figure 5b). The results showed that Acinetobacter, Enterococcus, Erysipelotrichaceae, and Ruminococcaceae_UCG_014 were the main microorganisms responsible for the differences between the two groups and were the main regulatory bacteria for 99Tc-MDP treatment.

Figure 5 LEfSe analysis. (a) The magnitude of the species’ contribution to the difference. LDA SCORE indicates the magnitude of the difference contribution; the larger the absolute value, the greater the difference contribution. (b) Cladogram annotation of the differential species. Nodes at each layer indicate phylum/class/order/family/genus from the inside out. Node diameter is proportional to relative abundance. Red, rich in the Before group; Green, rich in the After group; Yellow, the difference is not significance.

Random Forest Analysis to Identify Key Classified Bacteria

Random Forest model can efficiently and accurately classify microbial community samples. Mean Decrease Gini is used to assess the importance of variables. A higher value indicates that the variable (species) is more important. The results showed that Bacteroides, Escherichia-Shigella, and Prevotella_9 were the most important variables in the classification model and were the key bacteria to distinguish 99Tc-MDP before and after treatment (Figure 6a).

Figure 6 Targeted bacterial analyses. (a) Random Forest analysis. Mean Decrease Gini assesses the importance of the species, with higher values indicating more important species. (b) Association analysis. Red, positive correlation; Blue, negative correlation; *P < 0.05; **P < 0.01; ***P < 0.001.

Association Analysis of Gut Microbiota

Association analysis was performed on the top 30 genera with relative abundance in the gut microbiota. As shown in Figure 6b, Bacteroides was negatively correlated with Prevotella_9 and Ruminococcaceae_UCG-002. Prevotella_9 was positively correlated with Alloprevotella and Prevotella_2.

Functional Enrichment Analysis

KEGG pathway enrichment analysis was performed at three levels. At level 1 (Figure 7a), microbial functions in Metabolism, Genetic Information Processing, Cellular Processes, Human Diseases, Environmental Information Processing, and Organismal Systems were significantly different (P < 0.05) after 99Tc-MDP treatment. At level 2 (Figure 7b), there were 41 significantly different functions between the two groups, of which 43% (18 items) were Metabolism and Genetic Information Processing related. At level 3 (Supplementary Table 1), we found 293 differential functional pathways of which 30% (88 items) were related to Metabolism. COG database was used to infer the function of microbiota. The results showed that Carbohydrate transport and metabolism, Transcription, Amino acid transport and metabolism were the functions with the largest percentage (Figure 7c).

Figure 7 Functional enrichment analysis. (a) KEGG functional pathway at level 1. (b) KEGG functional pathway at level 2. (c) COG Function Category.

Discussion

The human microbiome is involved in the development and regulation of the host immune system, and dysbiosis of the intestinal flora is considered a key predisposing factor for RA.^23,24 An increasing number of studies have showed significant differences in the composition of gut microbiota in patients with RA compared to healthy controls. In China, 99Tc-MDP is often combined with MTX for the treatment of RA with anti-inflammatory, analgesic and immunomodulatory effects.^17,25,26 However, the relationship between 99Tc-MDP and gut microbes has not been reported. Therefore, we clarified the effect of 99Tc-MDP on the gut microbiota by 16S rRNA gene sequencing.

Consistent with recent reports, Bacteroidetes, Firmicutes, and Proteobacteria are the main dominant phylum in patients with RA and show an increase in Bacteroidetes and a decrease in Firmicutes.^7,27,28 A study on the prediction of MTX efficacy showed that the ratio of Firmicutes/Bacteroidetes in MTX non-responder is high.¹⁷ However, after 99Tc-MDP treatment, Bacteroidetes increased, Firmicutes decreased, and the ratio decreased. This suggests that 99Tc-MDP and MTX may have a synergistic effect in regulating intestinal flora and that 99Tc-MDP treatment can increase the response rate of MTX. Notably, our study found that 99Tc-MDP treatment did not alter gut microbial diversity in RA patients. This suggests that 99Tc-MDP is not simply an “antibiotic” or a “culture medium”, but specifically regulates the target flora to alleviate RA.

In patients with RA, there are large amounts of secreted immunoglobulins (IgA and IgM) such as RF, ACPA, and anti-carbamylated protein antibodies. Studies have found that these immunoglobulins are usually produced in mucosal areas and are associated with Prevotella.²⁹ An HLA-DR presenting peptide was identified in Prevotella’s 27-kDa protein (Pc-p27) that stimulated Th1 cell responses in 42% of patients with new-onset RA.³⁰ It was also shown patients with RA (both new-onset and chronic) either had IgA-like antibody responses to Pc-p27 or Prevotella or presented anti-Prevotella IgG antibodies.⁴ Random Forest analysis also showed that Prevotella is a key classified bacterium.

Enterococcus can have a broad inhibitory effect on pathogenic bacteria by producing bacteriocins.³¹ Enterococcus may be a potential probiotic for RA patients, enhancing the anti-inflammatory effects of MTX and alleviating arthritic symptoms in adjuvant-induced arthritis rat model.³² Our study showed a significant increase in Enterococcus after 99Tc-MDP treatment, while LEfSe analysis also identified Enterococcus as the main differential bacteria. Therefore, the regulation of Enterococcus abundance may be one of the mechanisms of 99Tc-MDP for RA treatment. In addition, several studies have found an increased proportion of Escherichia and Ruminococcaceae in fecal samples from RA patients.^33,34 And 99Tc-MDP treatment can reverse the regulation and reduce the abundance of these bacteria.

Functional enrichment analysis showed the effects after 99Tc-MDP treatment were focused on the Metabolism and Genetic Information Processing functions of gut microbiota. Short-chain fatty acids (SCFAs) are metabolites of intestinal flora. While providing energy for the host, SCFAs also promote the activation of Tregs by inhibiting histone deacetylase (HDAC).^35,36 According to studies, levels of valproic acid, butyric acid and propionic acid were reduced in RA.^37–39 However, a high-fibre diet increases these SCFAs, reduces pro-inflammatory cytokines (IL-33, IL-18 and CCL2), and prevents inflammation-induced bone loss. This coincides with the role of 99Tc-MDP in bone metabolism. In addition, overactivation of tryptophan metabolism contributes to RA, and dysregulation of tryptophan metabolism is closely associated with Escherichia.^33,40–43 On the other hand, the gut microbiome also interacts with many RA-related genes. RA susceptibility gene HLA-DRB1 favors a more inflammatory gut microbiome that drives aberrant immunity in RA.^6,44,45 The relative abundance of Prevotella in RA is negatively correlated with the presence of HLA-DRB1.⁴⁶

Conclusions

In conclusion, this study revealed differences in the composition of gut microbiota before and after 99Tc-MDP treatment. 99Tc-MDP exerts therapeutic effects by modulating the abundance of specific bacteria, especially Bacteroidetes, Prevotella, and Enterococcus. Influencing metabolism and genetic information processing of gut flora may be one of the mechanisms of how 99Tc-MDP treats RA.