Butyrate producing microbiota are reduced in chronic kidney diseases

Chronic kidney disease is a major public health concern that affects millions of people globally. Alterations in gut microbiota composition have been observed in patients with chronic kidney disease. Nevertheless, the correlation between the gut microbiota and disease severity has not been investigated. In this study, we performed shot-gun metagenomics sequencing and identified several taxonomic and functional signatures associated with disease severity in patients with chronic kidney disease. We noted that 19 microbial genera were significantly associated with the severity of chronic kidney disease. The butyrate-producing bacteria were reduced in patients with advanced stages of chronic kidney diseases. In addition, functional metagenomics showed that two-component systems, metabolic activity and regulation of co-factor were significantly associated with the disease severity. Our study provides valuable information for the development of microbiota-oriented therapeutic strategies for chronic kidney disease.

patient groups, meanwhile the CKD5 group showed a significant decrease in Simpson diversity compared with CKD 3A (FDR < 0.06), CKD3B (FDR < 0.06) and CKD4 (FDR < 0.07) groups (Fig. 1A). The variation within patient groups was greater than variation across groups, and the overall gut microbiota composition of different CKD patient groups was not statistically different (Adonis test p-value = 0.057, Fig. 1B) even though ESRD patients tended to group separately from the patients with earlier stage patients (CKD 3A/3B, Fig. 1B). Mean Bray-Curtis distance within the cohorts was shown to be higher in the later stages of CKD (CDK 4/5 and ESRD) when compared with samples from subjects with early stage disease (CKD 3A/3B) (Fig. 1C). Mean Bray-Curtis distance of samples from the earliest stage group (CKD 3A) showed that the microbial composition of the subjects in the ESRD group is farthest away from patients in CKD 3A stage (Fig. 1D). Spearman correlation analysis showed that the dissimilarity in BMI between the subjects did not show any significant correlation with the taxonomic dissimilarity between the subjects (p-value = 0.182, Supplementary Fig. S1). In addition, the HIV status of the subjects did not have any significant effect on the taxonomic dissimilarity (ANOVA FDR > 0.19, Supplementary Fig. S2).

Gut microbiota composition.
A total of 24 microbial species were significantly different between early stage patients (CKD 3A/B) and ESRD patients (FDR < 0.1, Fig. 2A). A total of 57 microbial species were significantly associated with the disease severity (Supplementary Table S2). The top six positive and negative associations are shown in Fig. 2B, among which Eubacterium rectale and species belonging to Collinsealla genera demonstrated the strongest association. At genus level, 19 microbial genera were significantly associated with the severity of CKD (Supplementary Table S3), and the top six positive and negative correlations are shown in Fig. 2C.
Further, we identified 22 microbial species harboring genes encoding butyrate kinase and butyryl-coA transferase in patients. The butyrate synthesis associated species composition was estimated by summing the relative composition of 22 known butyrate synthesis species identified in the subjects at > 1% relative composition. The summed up relative composition from these species represents the butyrate synthesis capacity. The Linear mixed effect model showed the butyrate synthesis associated taxa was negatively associated with the disease severity (p-value = 6e −04 ) and ESRD patients had the lowest level of butyrate synthesis associated taxa (Fig. 3A). Spearman correlation analysis showed that the serum creatinine levels were significantly correlated with butyrate synthesis associated species composition (Fig. 3B, Left). The albumin/creatinine ratio in the serum showed significantly negative association with butyrate synthesis associated species composition (Fig. 3B, Right). The level of Bifidobacterium was increased with disease severity, meanwhile lactobacillus was decreased with disease severity (Fig. 3C). Interestingly, Methanobacteria increased from early to late stage CKD patients, but it was not detected in patients with ESRD (Fig. 3D).
Next, we identified five microbial genera and eight microbial species that were different between the diabetic patients and non-diabetic patients (FDR < 0.01, Fig. 4A,B). However, the diabetic status didn't significantly affect the gut microbial community structure (Adonis test p-value > 0.874, Supplementary Fig. S3). Eight microbial genera and 11 microbial species were different between patients with high (albumin/creatinine ratio > 300) and low (albumin/creatinine ratio < 30) albuminuria (FDR < 0.01, Fig. 4C,D).

Discussion
Using shot-gun metagenomics, we identified microbial taxa and functional signatures associated with the disease severity in patients with CKD. Notably, butyrate producers reduced with disease severity (Fig. 3A). Butyrate is a four-carbon short-chain fatty acid, produced by the microbial digestion of dietary fibers. Butyrate is a major energy source for colonocytes 10 . Butyrate interacts with G protein-coupled receptors (GPCRs) and acts as an inhibitor for histone deacetylase (HDAC). Butyrate has anti-inflammatory effects, enhances intestinal barrier function and mucosal immunity 11 . The decrease in butyrate producers in advanced stage of CKD might have adverse effects on the host. In addition, Methanobacteria increased from early to late stage CKD patients (Fig. 3D). We have demonstrated that intestinal colonization with methanogenic Archaea lowers the Trimethylamine N-oxide level in plasma, with a tendency to attenuate the atherosclerosis burden in Apoe −/− mice 12 . Several Collinsella species were positively associated with the CKD disease severity (Supplementary Table S2). Increased abundance of Collinsella has been reported to be associated with increased cholesterol levels in healthy   14 . Bacteria sense, respond and adapt to the changes and stimuli in their environment through two-component systems, which contain a sensor protein-histidine kinase for receiving external input signals and a response regulator that conveys a proper change in the bacterial cell physiology 15 . Among the 17 significant KEGG modules associated with disease severity, 13 were two-component regulatory systems (Fig. 5B). These two-component regulatory systems were involved in stress response (CpxA-CpxR, LiaS-LiaR), antimicrobial peptide resistance (BasS-BasR), capsule synthesis (RcsC-RcsD-RcsB), chemotaxis (CheA-CheYBV), signaling (RpfC-RpfG), cell fate control (PleC-PleD), metabolism (CitA-CitB, DcuS-DcuR, UhpB-UhpA), and regulation of phosphate, sodium and potassium (CreC-CreB, KdpD-KdpE, NatK-NatR). These changes in the two-component systems suggest that the gut microbiota respond and adapt to the changes in the environmental niches when early stage CKD progressed to ESRD.
The metabolic capacity of the gut microbiota also differed when CKD progressed to advanced stage. Four KEGG module involved in different compound metabolism were negatively correlated with the disease severity, including hexose phosphate uptake, trehalose biosynthesis, citrate fermentation and catechol meta-cleavage, one of the major pathways for the degradation of aromatic compounds. In addition, among five KEGG enzymes negatively correlated with the disease severity, three were involved in the metabolism, including 3-(hydroxyamino) phenol mutase, glucose-1-phosphatase and dimethylmaleate hydratase large subunit. 3-(hydroxyamino)phenol mutase is involved in the degradative pathway of 3-nitrophenol, a phenolic compound. Glucose-1-phosphatase primarily acts as a scavenger for glucose 16 . Dimethylmaleate hydratase catalyzes the conversion of (2R,3S)-2,3-dimethylmalate to dimethylmaleate 17 , which is a dimethyl ester of maleic acid.
In addition to the metabolic activity, the gut microbiota's functional capacity of cofactor biosynthesis and regulation of antibiotics were also affected by the disease severity. Molybdopterin-synthase adenylyltransferase and streptothricin hydrolase were negatively correlated with the disease severity (Fig. 6B). Molybdopterinsynthase adenylyltransferase participates in the synthesis of molybdopterins, which are a class of cofactors found in a large group of enzymes. Streptothricin hydrolase catalyzes the hydrolysis of streptolidine lactam group in streptothricin antibiotics and inactivate them 18 . Capreomycidine synthase was positively associated with the disease severity (Fig. 6B), which catalyzes the biosynthesis of capreomycidine, a nonproteinogenic amino acid which is involved in the biosynthesis of the tuberactinomycin class of peptide antibiotics such as viomycin and capreomycin.
There are some limitations of this study. The sample size is relatively small, and the study population is heterogeneous with relatively small numbers of patients in each CKD stage. Thus, the findings from this study need to be validated in a larger cohort. Treatment information was not available in this patient cohort, which could influence the intestinal microbiota. Although we found the metabolic capacity of the gut microbiota altered with the disease severity, its impact on the host metabolism needs to be further validated using metabolomics approaches. In addition, further experiments are needed to investigate the causal relationship between the change of gut microbiota and CKD, such as fecal microbiota transplantations. In conclusion, this study provides valuable information about the association between disease severity and the dysregulation of the gut microbiota composition and its functional capacity. This knowledge can be valuable in the development of gut microbiota targeted precision medicine for patients with different stages of CKD.  DNA sequencing of fecal bacterial cell pellets and microbiome analysis. The DNA of fecal bacterial pellets were extracted and sequenced at Diversigen using BoosterShot Shallow Shotgun Sequencing. A mean sequencing depth of ≥ 0.5 M reads were used per sample. Samples were extracted with MO Bio PowerFecal or MO Bio PowerSoil Pro (Qiagen) automated for high throughput on QiaCube (Qiagen), with bead beating in 0.1 mm glass bead plates. Samples were then quantified using Quant-iT Picogreen dsDNA Assay (Invitrogen). a procedure adapted from the Nextera Library Prep kit (Illumina) was then used for preparing the libraries. Libraries were sequenced either on an Illumina NextSeq using single-end 1 × 150 reads with a NextSeq 500/550 High Output v2 kit (Illumina). DNA sequences were first filtered for low quality (Q-Score < 30) and length (< 50) followed by trimming of adapter sequences using cutadapt 20 software. Fastq files were then compiled into a single fasta using shi7. Sequences were then trimmed to a maximum length of 100 bp. The filtered and trimmed DNA sequences were then aligned to a database containing all representative genomes in RefSeq for bacteria with additional manually curated strains compiled by Diversigen. The reads were aligned at 97% identity against the reference genomes in the database. All sequences were compared to all reference sequence in the curated database using fully gapped alignment with BURST 21 . Ties were broken by minimizing the overall number of unique Operational Taxonomic Units (OTUs). Each input sequence was assigned the lowest common ancestor that was consistent across at least 80% of all reference sequences tied for best hit for taxonomy assignment. The OTU tables were rarified at 10,000 reads and samples with fewer than 10,000 reads were discarded. KEGG functional profiling was obtained by using the method similar to one previously described 22 . The profiles were summarized at KEGG Ortholog (KO) levels and KEGG Module levels.
Stool butyrate measurement. Due to the sample availability, butyrate concentrations were measured in stool samples collected from 8 patients with stage 3A, 6 patients with stage 3B, and 8 patients with stage 4 and 5. Stool butyrate measurement was performed at the West Coast Metabolomics Center following standard protocol.

Statistical analysis.
To evaluate effect on the overall microbiome structure, we square-root transformed relative abundances of bacterial species and calculated pair-wise bray-curtis dissimilarities between the subjects. We then performed principal coordinate analysis on the dissimilarity matrix. Further, we performed Wilcoxon rank test for identifying significantly different taxa between (1) early stage CKD (3A, 3B) and ESRD patients, (2) diabetic and non-diabetic patients and (3) high and low albuminuria patients. The test results for all comparisons were corrected using FDR. Only taxa with relative abundances greater than or equal to 0.1% in either group were examined. The same approach was used to identify significantly different KEGG modules and KO families between (1) early stage CKD (3A, 3B) and ESRD patients. We used a heuristic scoring scheme for representing the severity of the subjects. CKD 3A subjects were assigned a severity score of 10, CKD 3B subjects were assigned 20, CKD 4 subjects were assigned 30, CKD 5 stage subjects were assigned 40 and ESRD subjects were assigned a score of 50. MaAsLin2 was then used to explore the effect of disease severity on the microbial relative abundance and functional composition of the microbial metagenome 23 . FDR < 0.1 was considered as significant.