We performed a cross-sectional collection of fecal samples from 23 AAV patients with kidney injury, 15 patients with lupus nephritis (LN), and 27 healthy controls in West China Hospital of Sichuan University. Participants with AAV were eligible if they met the modified American College of Rheumatology classification criteria (Jennette et al., 1994; Seeliger et al., 2017). Patients with LN fulfilled the Systemic Lupus International Collaborating Clinics (SLICC) classification criteria and some proven by renal biopsy (Hochberg, 1997; Weening et al., 2004). Fecal samples from 27 healthy donors were collected as healthy controls. Disease activity was evaluated using the clinical composite Birmingham Vasculitis Activity Score (BVAS) V3.0 for patients with AAV and Systemic Lupus Erythematosus Disease Activity Index (SLE-DAI) for those with LN (Bombardier et al., 1992; Suppiah et al., 2011). All participants with AAV, LN, and healthy controls were excluded if they were younger than 18 years or had another systemic inflammatory disorder, known history of HIV or primary immunodeficiency, lymphoma, or other diseases that mimics AAV. Fresh fecal samples were collected in sterile containers and immediately stored at liquid nitrogen within a maximum of three hours from defecation. They were subsequently transported to the laboratory and stored at −80°C until DNA extraction. Written informed consent was obtained from all the participants in the study. Biomedical Ethics Committee of West China Hospital of Sichuan University approved the study (No. 2016–273).

Microbial DNA was extracted from 65 fecal samples using the E.Z.N.A.^® soil DNA Kit (Omega Bio-tek, Norcross, GA, U.S.) according to the manufacturer’s protocols. The V3-V4 regions of the bacterial 16S rRNA gene sequences were amplified from the diluted DNA extracts with the primers 338F (5′- ACT​CCT​ACG​GGA​GGC​AGC​AG-3′) and 806R (5′-GGACTACHVGGGTWTCTAA T-3′). Polymerase Chain Reaction (PCR) amplification was performed in a 20 μL mixture containing 4 μL of 5 × FastPfu Buffer, 2 μL of 2.5 mM dNTPs, 0.8 μL of each primer (5 μM), 0.4 μL of FastPfu Polymerase and 10 ng of template DNA. The reactions were hot-started at 95°C for 3 min, followed by 27 cycles of 95°C for 30 s, 55 °C for 30 s, and 72°C for 45 s, with a final extension step at 72°C for 10 min. The resulted PCR products were extracted by a 2% agarose gel and further purified using the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA) and quantified using QuantiFluor™-ST (Promega, USA) according to the manufacturer’s protocol.

Purified amplicons were conducted using Illumina MiSeq platforms (Illumina Inc., San Diego, CA) according to the standard protocols by Majorbio Bio-Pharm Technology Co. Ltd. (Shanghai, China). Before assembly, sequences were firstly filtered to remove low-quality or ambiguous reads, including reads with an average quality score <20 over a 50 bp sliding window, reads lacking exact matching with the primer, and reads containing ambiguous character. Operational taxonomic units (OTUs) were clustered with 97% similarity cut off using UPARSE (version 7.1 http://drive5.com/uparse/) and chimeric sequences were identified and removed using UCHIME. The taxonomy of each 16S rRNA gene sequence was analyzed by the RDP Classifier algorithm (http://rdp.cme.msu.edu/) against the Silva (SSU123) 16S rRNA database using a confidence threshold of 70%.

To test whether gut microbial species could be differentiated among AAV patients with kidney injury, patients with LN, and healthy controls, a metric multidimensional scaling method based on projection known as principal coordinates analysis (PCoA) and the Adonis test were used. We also compared alpha diversity indices across these three groups, such as Sobs, Chao, Shannon, and Simpson., Linear discriminant analysis Effect Size (LEfSe) analysis was used to identify differentially abundant bacteria among the three groups with a cutoff of 2.0. Microbial phylum and genera concentrations among AAV patients with kidney injury, patients with LN, and healthy controls were compared using a Mann-Whitney test. Using a normalized relative abundance matrix, LEfSe performs the Kruskal-Wallis rank sum test to determine the features with significantly different abundances between assigned taxa and uses LDA to assess the effect size of each feature (Segata et al., 2011). Spearman’s correlation was used to investigate the association between the levels of clinical indicators and fecal microbiota composition in AAV patients with kidney injury and visualized by heatmap in R using the “heatmap” package. The PICRUSt was used to construct the Kyoto Encyclopedia of Genes and Genomes (KEGG) Orthology (KO) and KEGG pathway/module profile, predicting functional profiling of microbial communities using 16S rRNA marker gene sequences.

We continuously enrolled 23 AAV patients with kidney injury, including 15 females and 8 males with a median age of 62 years in this study. Among the 23 AAV patients, two were PR3-ANCA positive, twenty were MPO-ANCA positive, and one was both PR3-ANCA and MPO-ANCA positive. The median estimated glomerular filtration rate (eGFR) and serum creatinine were 12.5 ml/min/1.73 m² and 349 μmol/L, respectively. The clinical characteristics of the AAV patients with kidney injury, the patients with LN, and the healthy controls were summarized in Table 1.

As shown in Figure 1A, the alpha-diversity indexes, including Sobs, Shannon, and Chao indexes, were significantly lower in AAV patients with kidney injury than those of healthy controls (Sobs P < 0.001, Shannon P < 0.001, Chao P < 0.001). At the same time, there were no significant differences in the above indexes between AAV patients with kidney injury and patients with LN. The beta-diversity difference using the Bray Curtis was exhibited in PCoA analysis (Figure 1C), which demonstrated a clear difference among AAV patients with kidney injury, patients with LN, and health controls (ANOSIM, p = 0.001).

As exhibited in the Venn diagram (Figure 1B), the number of Operational Taxonomic Units (OTUs) in common between healthy controls and AAV patients with kidney injury was 636, while AAV patients with kidney injury and patients with LN had 524 OTUs in common. The group of AAV patients with kidney injury had 58 specific OTUs, and the group of patients with LN had 12 specific OTUs, which the other two groups did not share. To determine different taxa from phylum to genus level among the three groups, the linear discriminate analysis (LDA) effect size (LEfSe) algorithm was also employed. The LEfSe showed two classes, two orders, one family, one genus were higher in feces of AAV patients with kidney injury and one family, two genera were higher in feces of patients with LN, in comparison with that of healthy controls (Figures 1D,E).

Taxonomic assignment of the OTUs helped to reveal the composition of bacterial population from phylum to genus level separately. At the phylum level, Firmicutes was the most predominant bacteria of gut microbiota, followed by Bacteroidetes and Proteobacteria in all groups. Firmicutes and Bacteroidetes had a lower proportion in feces of AAV patients with kidney injury and those of LN patients, compared to those of healthy controls. Whereas the abundance of Proteobacteria was higher in the feces of AAV patients with kidney injury and LN than that of healthy controls (Figure 2A). At the class level, significant differences were observed in the bacterial composition; the ascendant bacteria in feces of three groups were Clostridia, Bacteroidia, Bacilli, Negativicutes, and Gammaproteobacteria (Figure 2B). The proportions of Clostridia, Negativicutes, Betaproteobacteria, and Fusobacteriia were lower, while Bacilli and Verrucomicrobiae were higher in samples of AAV patients with kidney injury compared with those of healthy controls. The decreased species assigned to Negativicutes and Betaproteobacteria as well as the increased species assigned to Bacilli in LN group were compared with healthy controls (Figure 2D). Notably, prominent differences in gut microbiome composition were observed at the family level. The abundance of the bacterial family Enterobacteriaceae, Streptococcaceae, and Verrucomicrobiaceae were higher in samples from AAV patients with kidney injury than healthy controls. Conversely, the relative abundance of family Lachnospiraceae, Prevotellaceae, and Ruminococcaceae were lower in samples from AAV patients with kidney injury compared with those from healthy controls. (Figures 2C,E). At the genus level, the abundance of Enterococcus, Peptoclostridium, Streptococcus, and Akkermansia was higher in AAV patients with kidney injury compared with healthy controls. However, the proportion of Faecalibacterium, Prevotella_9, and Roseburia, was lower in samples of AAV patients with kidney injury compared with those of healthy controls (Table 2).

n order to evaluate whether the proportion of changing gut microbiota was related to serum creatinine (SCR), blood urea nitrogen (BUN), eGFR, and hemoglobin (Hb), the Spearman’s correlation analysis was employed (Figures 3A–C). Among the genera analyzed, at the class, family, and genus level, a negative association between the level of SCR and the proportion of Deltaproteobacteria, unclassified_o_Bacteroidales, Prevotellaceae, Desulfovibrionaceae Paraprevotella, and Lachnospiraceae_NK4A136_group was observed. In contrast, the proportion of Deltaproteobacteria, unclassified_o_Bacteroidales, Desulfovibrionaceae, Paraprevotella, and Lachnospiraceae_NK4A136_group had positive correlations with the level of eGFR.

Compared with LN patients, based on Bray-Curtis dissimilarity at OTU level showed that the microbiota composition in samples of AAV patients with kidney injury was similar to that of LN patients (ANOSIM, p = 0.653, Figure 4A). In addition, the difference of gut microbiota composition between AAV patients and LN patients was explored using the Mann-Whitney test at different taxon level. There was no apparent difference in the phylum, class, and family levels between samples from AAV patients and LN participants (Supplementary Figures 1A–C). Instead, the proportion of order Pseudomonadales was significantly higher in samples of AAV patients with kidney injury patients than that of patients with LN (Figure 4B).

By variance analysis of KEGG metabolic pathways, the differences of metabolic pathways of functional genes in microbiota between the samples of AAN patients with kidney injury and healthy controls could be observed (Figure 4C). Between these two groups in the second-level KEGG pathways, many pathways like “Carbohydrate metabolism”, “Signal transduction”, “Lipid metabolism”, and “Immune disease” were lower in AAV patients with kidney injury than that of healthy controls.