Changes of the gut microbiota composition and short chain fatty acid in patients with atrial fibrillation

Study population

This study enrolled 60 patients who hospitalized in the department of cardiology the First Affiliated Hospital of Wenzhou Medical University from January 2020 to December 2020, including 30 patients in the sinus rhythm (SR) group and 30 patients in the AF group. The study protocol was approved by the Ethics Committee in Clinical Research of the First Affiliated Hospital of Wenzhou Medical University, and all study participants provided written informed consent. The approval reference number was clinical research ethics approval number #174. Portions of this text were previously published as part of a thesis (Qu et al., 2019; Qu et al., 2021).

Patients with SR and AF were grouped as previously described by Qu et al. (2019). Specifically, patients were included in the SR group if their electrocardiogram (ECG) suggested sinus rhythm. In addition, SR patients had not received any previous antiarrhythmic treatment, had no history of AF, and had no evidence of AF in the results of electrocardiogram, dynamic electrocardiogram, or ECG telemetry during hospitalization. Patients were categorized to AF based on previous medical history and ECG, dynamic ECG or cardiac telemetry system results during hospitalization. AF was defined as cardiac arrhythmia with “absolutely” irregular RR intervals and no distinct P waves on ECG. The main exclusion criteria were patients with other types of arrhythmias eventually determined by ECG or electrophysiological examination; acute and chronic heart failure patients with New York Heart Association level III–IV requiring hemodialysis or ventilator therapy; liver function incomplete with transaminase was more than upper limit three times; acute or chronic renal insufficiency with creatinine more than three times of the upper limit; patients with acute myocardial infarction with Killip III-IV requiring ventilator support or intra-aortic balloon dilation; patients with hemodynamic instability due to various reasons requiring rescue; patients within one month prior to the enrolled had intestinal diseases or those who had taken probiotics or antibiotics, because these factors may affect GM and SCFA in faeces; patients did not consent to participate in this study or any other reasons that the investigator felt that the patient’s condition would interfere the clinical study. SR and AF’s ECG were blindly adjudicated by a validation committee who referred to the Guidelines for the Management of Atrial Fibrillation (ESC) (Kirchhof et al., 2016).

Related baseline and clinical samples

The data were collected from the hospital database using a web-based electronic medical record system. Data were collected as previously described by Qu et al. (2021). Specifically, the following demographic and clinical characteristics were extracted: age, gender and body mass index (BMI, kg/m²). Traditional cardiovascular risk factors including smoking, alcohol intake, hypertension, diabetes and coronary heart disease. key laboratory data considered as prognostic markers were also collected. Including C-reactive protein (CRP), eGFR (calculated based on serum creatinine), alanine aminotransferase (ALT), aspartate aminotransferase (AST), low-density lipoprotein cholesterol (LDL-C), brain natriuretic peptide (BNP). Baseline peripheral vein blood samples were collected on the morning following hospitalization from patients in stable condition. All the biochemical measurements were performed at the hospital’s laboratory.

Stool samples collection

Stool samples were collected from middle section of the fresh feces for the first time after hospitalization. Stool samples were divided into well-sealed 250 mg sterilized tubes by the patients themselves, according to our instructions, and stored the samples at −20 °C. Researchers collected the stool samples within 24 h, frozen with liquid nitrogen for 15 min, then stored at −80 °C until analysis.

DNA extraction and PCR amplification

DNA was extracted from the stool samples and 16S rRNA amplicon sequencing were performed using the NovaSeq PE250 sequencing platform, according to the manufacturer’s instructions. Data were collected as previously described in Chen et al. (2020). Specifically, Genomic DNA of the samples were extracted by SDS-CTAB method and the purity and concentration of DNA were detected by agarose gel electrophoresis. Appropriate amount of DNA was taken and diluted to 1ng/µl in sterile water. Using diluted genomic DNA as a template, the V4 region of the bacterial 16S rRNA was amplified by PCR with the specific primers 515F(5′-GTGCCAGCMGCCGCGGTAA-3′) and 806R(5′-GGACTACHVGGGTWTCTAAT-3′) labeled in a 12 bp barcode. Using high-fidelity PCR Master Mix with GC Buffer (Phusion^®; New England Biolabs, Ipswich, MA, USA) and high-fidelity enzyme to ensure amplification efficiency and accuracy.

Sequencing and data processing

Reads for each sample were demultiplex from the raw data according to Barcode sequence and PCR primer sequence. Original tags data (raw tags) were formed by splicing the reads from each sample using FLASH (V1.2.7, http://ccb.jhu.edu/software/FLASH/) (Magoč & Salzberg, 2011). The spliced raw tags were strict filtered (Bokulich et al., 2013) to get high-quality clean tags according to published protocols (Qiime V1.9.1, http://qiime.org/scripts/split_libraries_fastq.html) (Caporaso et al., 2010). Clean tags detected chimeric sequences by comparing with the species annotation database (https://github.com/torognes/vsearch/) (Rognes et al., 2016), and finally removed the chimeric sequences to obtain the final effective tags. Operational Taxonomic Units (OTUs) were clustered against the Uparse software (Uparse V7.0.1001, http://www.drive5.com/uparse/) (Haas et al., 2011) at 97% identity for all effective tags. The representative sequence of the OTUs was analyzed by mothur method and SILVA138 (http://www.arb-silva.de/) (Edgar, 2013) SSU rRNA database (Wang et al., 2007). MUSCLE program (Version 3.8.31, http://www.drive5.com/muscle/) (Quast et al., 2013) was used for fast multi-sequence alignment to obtain the phylogenetic relationships of all OTUs. Finally, the data of all samples were homogenized, and the sample with the least amount of data was taken as the standard for homogenization. The subsequent Alpha and Beta diversity analysis were based on the homogenized data.

Alpha diversity analysis

We analyzed within-community microbial diversity by using alpha diversity analysis. Qiime v1.9.1 was used to calculate Observed-species, Chao1, ACE, Shannon, Simpson, PD_whole_tree indexes. R version 2.15.3 (R Core Team, 2013) was used to plotted the rarefaction curve, rank abundance curve and species accumulation boxplot. The difference analysis between AF group and SR group of alpha diversity indexes were conducted by t-test and Wilcox test.

Beta diversity analysis

We analyzed the dissimilarities in the microbiomes by using beta diversity analysis. Qiime v1.9.1 was used to calculate unweighted unique fraction metric (Unifrac) between AF group and SR group. Principal co-ordinates analysis (PCoA) and non-metric multi-dimensional Scaling (NMDS) diagrams were drawn using R v2.15.3. We performed PCoA by using WGCNA, STATS and GGplot2 packages of R and performed NMDS analysis by using vegan package of R. Analysis of similarities (Anosim), multi-response permutation procedures (MRPP) and Adonis (permutational multivariate analysis of variance) were analyzed by using the Anosim function, MRPP function and Adonis function of R vegan package respectively. The difference analysis in the unweighted Unifrac between the AF and SR groups were conducted by t-test and Wilcox test in R software.

LEfSe and PICRUSt analysis

LEfSe software was used for the linear discriminant analysis (LDA) effect size (LEfSe) analysis and the threshold of LDA score (log10) was set to 4. Permutation tests were performed between groups at various classification levels (Phylum, Class, Order, Family, Genus and Species) by using the metastats analysis of R software, and P-values were obtained. Then, the Benjamini and Hochberg false discovery rate method was used to correct P-values and obtained q-values (Edgar, 2004). The difference analysis between the AF and SR groups were also using the t-test and Wilcox test in R software. GM functional was predicted by using the phylogenetic investigation of communities by reconstruction of unobserved states (PICRUSt) (Langille et al., 2013). The detailed prediction process was referred to http://picrust.github.io/picrust/tutorials. Function prediction based on kyoto encyclopedia of genes and genomes (KEGG) database was performed based on 16S sequencing data. Additionally, Spearman correlation was used to associate abundant differential taxa with SCFA.

SCFA sample processing and gas chromatography analysis

The GCMS system mainly consisted of gas mass spectrometer (Thermo TRACE 1310-ISQ LT, USA), vortex mixer (QL-866, China) and refrigerated centrifuge (Xiangyi, China). The SCFA were extracted as described previously (Hoving et al., 2018). Briefly, 50 mg of the sample was added with 50 µL15% phosphoric acid, 100 µL of 125 µg/mL internal standard (isohexanoic acid) solution and 400 µL of diethyl ether homogenate for 1 min. The sample was centrifuged at 12,000 RPM for 10 min at 4 °C, and the supernatant was transferred to a chromatographic vial for GC analysis. Chromatographic conditions: agilent HP-INNOWAX capillary column (30m* 0.25 mm ID* 0.25 µm). The injection volume was 1 µl, and the split ratio was 10:1. The carrier gas was helium with a flow rate of 1.0 mL/min.

Statistical analysis

The statistics for continuous variables were expressed as means and standard deviations or as median and interquartile range when the distribution was not normal. Categorical variables as frequencies and proportions. Baseline characteristics were compared between AF group and SR group. A Fisher’s exact test or chi-square test was used to identify differences in categorical data, and independent-sample t test (data with normal distribution) or Mann–Whitney U test (data with non-normal distribution) were used to identify differences in quantitative data. Spearman correlation analysis was used to identify the correlation between fecal SCFA and GM. Positive R-values represented positive correlation and negative R-values represented negative correlation. Statistical analyses were performed using SPSS software version 26 (SPSS, Inc., Chicago, IL, USA). P values <0.05 were considered statistically significant.

The raw sequencing reads of this study have been deposited in National Population Health Data Center of China: https://www.ncmi.cn//phda/dataDetails.do?id=CSTR:17970.11.Z061U.202204.366.V1.0.

Baseline characteristics of the study cohort

In the current study, 60 patients consisting of 30 SR patients and 30 AF patients were included in our study. The baseline clinical characteristics of the two groups were shown in Table 1. The most baseline clinical characteristics among the SR and AF groups were similar. Compared to the SR group, AF patients were presented with higher BNP level.

Fecal SCFA levels in SR and AF groups

The concentrations of acetic acid, propionic acid, butyrate acid, isobutyric acid, valeric acid, isovaleric acid and caproic acid in feces are shown Table 2 and Fig. 1). The isobutyric acid level was significantly decreased in AF group (53.75 ± 37.49) compared with that in SR group (82.57 ± 52.23, P = 0.017). Besides, compared with SR patients (83.56 ± 70.66), the isovaleric acid was significantly decreased in AF patients (52.34 ± 43.50, P = 0.044). However, no significant differences were found in levels of acetic acid (P = 0.232), propionic acid (P = 0.653), butyric acid(P = 0.476), valeric acid (P = 0.655) and caproic acid (P = 0.511) between the two groups.

The differences of gut microbiota profiles between SR and AF groups

We drew the rarefaction curves to show that two groups approached a plateau, identifying that the sequencing depth was adequate in this study (Fig. S1A). Then we plotted the rank abundance curves to compare the richness and evenness between SR and AF groups and found reduced richness and evenness in AF patients (Fig. S1B). We further drew the species cumulation boxplot to determine the adequacy of sample size of our research and assess community richness (Fig. S1C). These curves indicated that majority of microbial diversity had been captured.

We also calculated the indexes between SR and AF patients to compare the richness by Sobs index, Chao1 and ACE (Figs. 2A, 2B, 2C), α-diversity by Shannon and Simpson (Figs. 2D, 2E), and phylogenetic diversity by PD whole tree (Fig. 2F). In the Sobs index, we showed unique OTU with specific relative abundance threshold >0.01% to filter out potential noise from sequencing error. These indexes suggested decreased species richness and α-diversity in AF patients compared to SR patients, but there were no statistical difference. However, the phylogenetic diversity by PD whole tree was significant decreased in AF group (P < 0.001). It indicated that SR patients had a higher abundance and diversity than in AF patients.

We then applied unweighted Unifrac analysis to compare similarities among GM communities (β-diversity). PCoA and NMDS diagrams revealed more clustered distribution from AF patients compared to SR patients. (Figs. 3A, 3B). Anosim analysis showed that the inter-group difference between AF group and SR group was significantly greater than the intra-group difference, indicating that grouping was meaningful (Fig. 3C. R = 0.098, P = 0.001). Both analysis suggested reduced β-diversity in AF group, which was further supported by MRPP analysis (A = 0.010, observed △ = 0.700,expected △ = 0.707, P = 0.001) and Adonis analysis (R² = 0.036, P = 0.001).

To show the difference in microbiome distribution between AF and SR cohorts for more direct visualization, we made the bar plot to showing the top 10 GM relative abundance of different patients by genus level (Fig. S2). To investigate the differential abundance of individual taxa, we first ran t-tests of the relative abundance of each taxon between SR and AF patients at different phylogenetic ranks, and found significant difference in four phylum, five classes, four orders, six families, five genera, and three species (Fig. 4). To analyze the statistical differences in microbial communities between the SR group and the AF group, we also performed the LEfSe analysis, which coupled statistical significance with biological consistency and effect relevance. The histogram reflected the LDA scores greater than 10⁴ that were computed for the features at the OTU level (Fig. 5A). Upon analyzing the fecal samples, a total of 5 species of bacteria were abundant in the SR group and 11 species were abundant in the AF group. The genera that were enriched in the AF group were Haemophilus_parainfluenzae, Haemophilus, Pasteurellaceae, Pasteurellales, Lactobacillus, Lactobacillaceae, Fusobacterium, Fusobacteriota, Fusobacteriales, Fusobacteriaceae and Fusobacteriia. Several genera were found in significantly high abundances in the SR group. These included Actinobacteria, Bifidobacterium, Bifidobacteriales, Bifidobacteriaceae and Dialister. Cladograms of the taxa with LDA values > 4.0 were depicted (Fig. 5B). Cladograms were generated from the LEfSe analysis, which showed the most differentially abundant taxa enriched in microbiota with green for the SR group and red for the AF group. Each small circle at different classification levels represented a classification at that level, and the diameter of the small circle was proportional to the relative abundance. The AF group showed a significant decrease in the class Actinobacteria, order Bifidobacteriales and family Bifidobacteriaceae and a greater abundance of the class Fusobacteriia, order Fusobacteriales, order Pasteurellales, family Lactobacillaceae, family Fusobacteriaceae and family Pasteurellaceae when compared with the SR group.

Spearman correlation was used to associate the differentially abundant taxa with the fecal levels of SCFA at the genus level (Fig. 6). Haemophilus which significantly increased in AF was negative correlated with isovaleric acid (rs =−0.27, P = 0.039) and isobutyric acid (rs =−0.26, P = 0.043). Other GM that was significantly increased in the AF and sinus rhythm groups did not differ significantly in SCFA.

Figure 3 was more about inter-group differences and intra-group aggregation, Fig. 3A was diagrammed by PCoA, while Fig. 3B was diagrammed by NMDS. PCoA and NMDS are mainly used to reflect the relation of sample distance matrix. The difference between the two lies in that, PCoA is based on the analysis of similarity distance between samples, while NMDS focuses more on reflecting the ordering relation of values in the distance matrix and weakens the absolute difference degree of values. In Fig. 3B, the smaller the Stress value is, the more complete information in the high-dimensional is retained. Figure 3C was based on the rank of Bray-Curtis distance value to test the significance of inter-group differences, which was used to test whether the inter-group differences were significantly greater than the intra-group differences. Fig. 4 was obtained by a t-test using abundance to identify species with significant differences (P < 0.05). Figure 5A mainly showed the species with significant difference abundance between the two groups. The histogram of distribution of LDA values showed species with an LDA Score greater than 4, namely the species with statistical difference between groups. Figure 5B was a Cladogram. The circles radiating from inside to outside represent taxonomic levels from phylum to species, which could more directly showed the species with significant difference abundance between the two groups.

Prediction of functional gene content

To realize the changes of intestinal microbial function in patients with AF, we based on KEGG database used PICRUSt. PICRUSt revealed six biological metabolism pathways at Level 1 pathways including metabolism, genetic information processing, environmental information processing, cellular processes, human diseases and organismal systems. Among them, metabolism, genetic information processing and environmental information processing dominated, accounting for 44.42%–48.17%, 17.29%–20.77% and 13.48%–17.92%, respectively. But, the six pathways in the AF group were not statistically different from those in the SR group. In KEGG pathway hierarchy level 2, we discovered 41 sub-functions. The top 10 were membrane transport, carbohydrate metabolism, amino acid metabolism, replication and repair, energy metabolism, translation, poorly characterized, metabolism of cofactors and vitamins, cellular processes and signaling, nucleotide metabolism. Within the two predicted functional categories at hierarchy level 2 pathway including amino acid metabolism (AF group 9.35% ± 0.33% vs SR group 9.58% ± 0.30%; P = 0.008) and metabolism of cofactors and vitamins (AF group 4.14% ± 0.13% vs SR group 4.25% ± 0.19%; P = 0.017) significantly decreased in AF patients. In KEGG pathway hierarchy level 3 we discovered 305 gene pathways. Oxidative phosphorylation(AF group 1.07% ± 0.06% vs SR group 1.11% ± 0.07%; P = 0.042), glycine, serine and threonine metabolism (AF group 0.97% ± 0.06% vs SR group 0.94% ± 0.05%; P = 0.024), photosynthesis proteins(AF group 0.79% ± 0.03% vs SR group 0.81% ± 0.03%; P = 0.047) and novobiocin biosynthesis (AF group 0.41% ± 0.06% vs SR group 0.44% ± 0.06%; P = 0.033) were significantly increased in the SR group. Whereas, naphthalene degradation(AF group 0.13% ± 0.02% vs SR group 0.12% ± 0.02%; P = 0.027) was significantly increased in the AF group. In KEGG pathway hierarchy level K we discovered 6309 taxa. K00145(N-acetyl-gamma-glutamyl-phosphate reductase: AF group 0.12% ± 0.02% vs SR group 0.13% ± 0.01%; P = 0.045) and K06148 (ATP-binding cassette: AF group 0.12% ± 0.02% vs SR group 0.11% ± 0.02%; P = 0.026) were the highest abundance gene pathways of GM in SR and AF patients, respectively

Study limitations

There were several limitations in this study. The most importantly, this study was an observational study and could not determine the causal relationship between GM or SCFA and AF; therefore, a prospective cohort study is needed to determine whether changes in GM or SCFA are associated with the etiology of AF. Furthermore, the sample size of this study was small, thus, larger studies are needed to validate our results. In addition, we cannot exclude the influence of confounding factors, such as sleep, diet and drugs, on the GM or SCFA of AF patients.