Alterations of gut fungal microbiota in patients with rheumatoid arthritis

Ethics statement

This study protocol was approved by the Ethics Committees of Dalian Municipal Central Hospital, Dalian, China (YN2021-003-01). All the procedures were performed in accordance with the guidelines approved by the Ethics Committee of Dalian Medical University, China. After receiving a written description of the aim of this study, all subjects signed written informed consent prior to study enrollment.

Subject recruitment and sample collection

From January 2020 to February 2021, a total of 62 patients with RA were recruited from Dalian Municipal Central Hospital in Dalian, China, and 29 healthy individuals also agreed to participate. RA patients fulfilled the RA criteria of the American College of Rheumatology (ACR)/European League Against Rheumatism (EULAR) 2010. Participants were excluded if they had taken antibiotics, probiotics or prebiotics within the previous 3 months, or had a history of inflammatory bowel disease or other autoimmune diseases, or gastrointestinal tract disorders. Fecal samples obtained from 42 RA patients and 39 healthy controls (HCs) were transported in liquid nitrogen and stored at −80 °C until fungal DNA was extracted. All serum samples were stored at −80 °C prior to use. Serum CRP, ESR, RF, IgG, IgA, IgM, IgE, CEA, CA19-9, and CA125 were measured on the day of stool sample collection in the laboratory at Dalian Municipal Central Hospital. An immunoturbidimetric assay was used to determine the values of serum CRP, RF, IgG, IgA, IgM, and IgE (Beckman Coulter, Brea, CA, USA). Serum CEA, CA19-9, and CA125 concentrations were determined by electrochemiluminescence (ECL) (Roche; Elecsys, Lenexa, KS, USA), and ESR was determined using an automatic erythrocyte sedimentation analyzer.

Genomic DNA extraction, PCR amplification and ITS2 gene sequencing

In accordance with the manufacturer protocol previously described, total metagenomic DNA in 200 mg (fresh weight) of each fecal sample was extracted using the QIAamp DNA stool minikit (Qiagen, Hilden, Germany) (Li, Leonardi & Iliev, 2019). The concentration and quality of DNA were measured by Nanodrop2000 (Thermo Fisher Scientific, Wilmington, CA, USA) and Agar Gel electrophoresis, respectively. The fungal microbiome was analyzed through the sequencing of internal transcribed spacer region 2 (ITS2). An ITS2 fragment (approximately 320 bases) was amplified with primers ITS3F (5′-GCATCGATGAAGAACGCAGC-3′) and ITS4R (5′-TCCTCCGCTTATTGATATGC-3′). The thermal cycle conditions were: initial denaturation at 98 °C for 2 min, followed by 25 cycles, denaturation at 98 °C for 15 s, annealing at 55 °C for 30 s, extension at 72 °C for 30 s, and extension at 72 °C for 5 min. After the individual quantification step, amplicons were pooled in equal amounts, and pair-end 2 × 300 bp sequencing was performed using the Illlumina MiSeq platform with MiSeq Reagent Kit v3 at Shanghai Personal Biotechnology Co. Ltd. (Shanghai, China).

Sequence analysis

The Quantitative Insights into Microbial Ecology (QIIME, v1.8.0) pipeline was employed to process the sequencing data, as previously described (Caporaso et al., 2010). Raw sequencing reads with exact matches to the barcodes were assigned to respective samples and identified as valid sequences. The low-quality sequences were screened using the following criteria (Gill et al., 2006; Chen & Jiang, 2014): sequences with a length less than 150 bp, sequences with mean Phred score less than 20, sequences with fuzzy bases, and single nucleotide repeats with a length greater than 8 bp. After chimera detection, the remaining high-quality sequences were clustered into operational taxonomic units (OTUs) with 97% sequence homology by UCLUST (Edgar, 2010). A representative sequence was selected from each OTU using default parameters. OTU taxonomic classification and OTU table generation was conducted as previously described (Zheng et al., 2020).

Bioinformatic analysis

The sequence data analysis was mainly performed using QIIME and R packages (v3.2.0). The community richness and diversity (alpha diversity) analysis (Chao1 and Shannon) was performed using QIIME 2 with the same sequence depth. The beta diversity analysis, principal coordinate analysis (PCoA), weighted or unweighted UniFrac analysis by R software were used to investigate the structural variation of microbial communities across samples. The linear discriminant analysis (LDA) effect size analysis (LEfSe) was performed to detect differentially abundant taxa across groups using the default parameters. The correlations between laboratory markers (including RF, CRP, IgG, IgA, IgM, IgE, CEA, CA19-9, CA125 and ESR) and fungal genera in fecal samples from RA patients were determined using Spearman’s method. A partial least squares discriminant analysis (PLS-DA) with a variable importance in projection (VIP) plot (Yuan et al., 2018) was performed to determine possible differences in OTUs between groups, which would help predict the functional contents of the metagenome. Key genera of VIP > 1.8 were considered important contributors to the model.

Statistical analysis

All data were presented as arithmetic mean ± standard error of mean (SEM). With the assistance of GraphPad Prism 7 (Graph Pad Software, La Jolla, CA, USA), the comparison between the two groups was statistically analyzed using the nonparametric t test, and multiple comparison adjustments were made with the FDR algorithm (Benjamini et al., 2001). Spearman’s rank correlation test was used to perform a correlation analysis. The statistical analysis was performed using SPSS version 19.0 (SPSS Inc., Chicago, IL, USA). Results were considered to be statistically significant with P < 0.05.

Accession number

The sequence data from this study were deposited in the NCBI Sequence Read Archive with the accession number RPJNA735021.

Investigation of clinical characteristics of RA patients and healthy controls

A total of 101 individuals were included in the current analysis. The clinical characteristics of all subjects are described in Table 1. The 62 (15 male, 47 female) patients with RA ranged in age from 33 to 80 years, with a mean of 60.8 years. The 39 (7 male, 32 female) healthy controls ranged in age from 32 to 80 years, with a mean of 60.4 year. Given that RA is three times more common in women than men (Favalli et al., 2019), there were more women than men in this study, which was consistent with the previous study (Sun et al., 2019). RA patients had higher CRP, ESR, RF, IgA, IgE, and CA125 levels than healthy controls, but there were no significant differences in other general biological indexes, such as IgG, IgM, CEA, and CA19-9 levels between the two groups (Table 1). RA patients were classified into subgroups by age and gender (Table 2), and the results showed that the levels of CRP and CEA in patients over 60 years old appeared increased compared with those under 60 years of age, while the level of IgM was decreased significantly in the older population (Table 2). However, these indexes showed no significant difference between the healthy controls when stratified by age (Table 3). The levels of IgA and IgE in male patients were significantly higher than those of female patients (Table 2), while IgE levels in the male control group showed a downward trend, but the difference was not statistically significant (Table 3).

Comparison of mycobiota between RA patients and healthy controls

In the present study, we obtained a total of 4,457,120 high-quality filtered reads from the fecal mycobiota of 39 HCs (48,412 ± 14,456) and 42 RA patients (61,167 ± 19,134) for subsequent analysis, with an average of 54,789 ± 18,112 reads per sample. The results showed that there were 526 OTUs unique to the HC group and 481 OTUs unique to the RA group, with a total of 1,138 OTUs in the two groups (Fig. 1A). Compared with HCs, the community richness (assessed using Chao1) and alpha diversity (assessed using the Shannon index) of intestinal fungi in RA patients was relatively decreased, but there were no statistical differences between both indexes (Figs. S1A and S1B, all P > 0.05).

Based on the Bray-Curtis, Jaccard, and Unweighted UniFrac algorithms, we then evaluated and compared the beta diversity of fungal flora in RA patients and healthy controls. The RA and healthy control groups could not be divided into different clusters (Adonis test, P > 0.05; Fig. S1C). Taken together, the alpha- and beta-diversity analyses demonstrated no significant changes in fungal community structure in RA patients compared with healthy controls. The analysis of the relative levles of fungi (Fig. 1B) showed that the Ascomycota and Basidiomycota phyla were the dominant taxa in both the HC and RA groups. However, there was no obvious change in the proportions of these two phyla in RA patients. At the genus level, classified genera including Aspergillus, Candida, and Penicillium were abundant in the fungal microbiota of both RA patients and healthy controls. Among the most dominant genera, Candida was increased in RA patients, but this increase did not reach statistical significance (P = 0.07, Fig. 1C). We found that the levels of Pholiota (P = 0.006), Scedosporium (P = 0.042), and Trichosporon (P = 0.037) in the RA group were lower than the HC group at the genus level (Fig. 1D).

LEfSe analysis was further performed to identify fungal groups with significant differences in levels between the RA and HC groups (LDA score > 2, P < 0.05; Fig. 1E). The comparison of the RA and HC groups revealed that most of the differential fungi could be classified at the genus level. Ganoderma, Auricularia, Pholiota, Aureobasidium, Pleurotus, Trichosporon, Podosphaera, Dipodascus, and Alternaria were enriched in the healthy controls, whereas Wallemia, Phanerochaete, Suhomyces, Tolypocladium, Trebouxia, Tremella, Scolecostigmina, and Sarocladium were prevalent in RA patients. Although most of these different genera were not abundant, our results suggested that fungal dysbiosis occurred in RA patients.

The Spearman correlation between mycobiota and physiological measures of RA was analyzed as well (Fig. 1F). The results indicated that the abundance of the vast majority of mycobiota prevalent in the healthy controls was negatively correlated with some serological features of RA; in contrast, the enrichment in RA patients was positively correlated with some serological characteristics of RA, except Wallemia, which was negatively related to the serum levels of IgG, IgA and CA125. As shown in Fig. S1D, the PLS-DA score plots showed model discrimination between RA patients and healthy controls. A total of three genera were recorded with VIP > 1.8 (Fig. S1E).

Comparison of mycobiota between RA patients aged less than 60 and over 60 years old

RA is an age-related disease, and RA patients were classified into two subgroups: those over 60 years old and those aged under 60. Sequencing results showed that there were 618 OTUs unique to the under 60 group and 280 OTUs unique to the over 60 group, with a total of 741 OTUs in the two groups (Fig. 2A). The calculated fungal richness and alpha-diversity indices showed no significant changes between the two groups (Figs. S2A and S2B, all P > 0.05). An analysis of beta diversity did not show a distinct clustering of samples in the under 60 and over 60 groups (Adonis test, P > 0.05; Fig. S2C). These data suggested no difference in the fungal composition of RA patients of different ages. The compositions of the fungal microbiota in RA patients under 60 years old and over 60 years old were assessed at the phylum and genus levels (Figs. 2B and 2C). We found that, at the genus level, the abundances of Dendroclathra (P = 0.009), Phacidium (P = 0.009), and Septobasidium (P = 0.009) in the Over 60 group were lower than those in the under 60 group (Fig. 2D). LEfSe analysis revealed Wallemia and Irpex as the most able to identify RA patients over 60 years old. The fungal microbiota that was able to distinguish RA patients under 60 years old were Olpidium, Pseudogymnoascus, Piromyces, Strobilomyces, and Septobasidium (Fig. 2E). Spearman correlation analysis demonstrated that Wallemia was negatively related to the serum level of IgG, Olpidium was negatively related to the serum level of CRP, Pseudogymnoascus was negatively related to the serum level of IgE, while Piromyces was positively related to the serum levels of IgA and IgM (Fig. 2F). PLS-DA score plot exhibited model discrimination between RA patients aged under 60 and over 60 years (Fig. S2D). A total of five genera were recorded with VIP > 1.8 (Fig. S2E).

Comparison of mycobiota between male and female RA patients

We classified RA patients into subgroups by gender and analyzed the sequencing results. In total, 694 OTUs overlapped between the male and female RA patients, in addition, 218 and 707 OTUs were uniquely present in the male and female RA patients, respectively (Fig. 3A). Analysis of alpha and beta diversities showed no significant changes between the two groups, suggesting no significant difference in fungal composition between female patients and male patients (Figs. S3A–S3C). The compositions of the fungal microbiota in male and female RA patients were analyzed at the phylum and genus levels (Figs. 3B and 3C). We found that the levels of Coryne (P = 0.006), Schizothecium (P = 0.006), Sarocladium (P = 0.006), Leotia (P = 0.006) and Scytinopogon (P = 0.044) were lower in female patients than in male patients at the genus level (Fig. 3D). LEfSe analysis revealed Pseudeurotiaceae was the most identifiable among female RA patients; while the fungal microbiota most able to identify male patients with RA were: Cladosporium, Coryne, Schizothecium, Arthrographis, Sarocladium, Eremomycetaceae, Raffaelea, Mycena, Leotia, Fusicolla, Scytinopogon, Mycoarthris, and Neurospora (Fig. 3E). Spearman correlation analysis demonstrated that Pseudeurotiaceae was negatively related to the serum level of IgE, Cladosporium was negatively related to the serum level of IgM, Arthrographis and Eremomycetaceae were negatively related to the serum level of CA19-9, while Fusicolla was positively related to the serum levels of RF, Sarocladium was positively related to the serum levels of IgE, and Neurospora was positively related to serum levels of IgA (Fig. 3F). A PLS-DA score plot revealed a clear discrimination between the microbial profile of female RA patients compared with male RA patients (Fig. S3D). Associated VIP scores (>1.8) allowed ranking of key microbial genera based on their importance in discriminating between female and male RA patients (Fig. S3E).