Study participants
The study protocols were approved by the Ethics Committee of Kaohsiung Medical University Hospital (KMUHIRB-E(I)-20160095 and KMUHIRB-E(I)-20180118) and Taipei Tzu Chi Hospital (07-X01-002). All participants provided written informed consent. Hemodialysis (HD) patients were recruited from the dialysis unit of Taipei Tzu Chi Hospital and Kaohsiung Medical University Hospital in Taiwan from August 2017 to February 2018. Eligible participants were those who received regular HD three times per week, 3.5–4 hours with high-flux dialyzers. Participants with active malignancies or participants who were prescribed antibiotics within three months before enrollment were excluded. Fecal samples were collected from 193 stable HD patients and analyzed by high-throughput 16S ribosomal RNA gene sequencing to compared participants with and without β-blocker treatment. All β-blocker users were prescribed for at least one month.
Comorbidity, laboratory, and clinical variables
All baseline characteristics of sociodemographic data, age, sex, dialysis vintage, arteriovenous shunt type, comorbidities, medications, and biochemical data were collected in the built-in electronic health care system. Blood samples were collected after overnight fasting through the arteriovenous fistula or graft before scheduled HD sessions. The biochemical data included serum values for hemoglobin, albumin, high sensitivity C reactive protein, total cholesterol, low-density lipoprotein, triglycerides, ion calcium, and phosphate from routine blood samples obtained within 30 days before stool sample collection. Diet was evaluated by a licensed dietitian using a modified short-form food frequency questionnaire. Diabetes was defined as HbA1C 6.5% or higher, or use of oral antidiabetic agents or insulin. Hypertension was defined as 140/90 mmHg or higher or taking blood pressure-lowering drugs. A history of myocardial infarction or documented by coronary angiography, class III or IV congestive heart failure, or a cerebrovascular accident were defined as cardiovascular disease.
Fecal sample collection and bacterial 16S rRNA amplicon sequencing and processing
All stool samples were frozen immediately after collection by each participant, then delivered in cooler bags to the laboratory (Germark Biotechnology, Taichung, Taiwan) within 24 hours. A QIAamp DNA Stool Mini Kit (Qiagen, MD, USA) was used to extract DNA from fecal samples. Barcode-indexed PCR primers (341F and 805R) were used to create an amplicon library by amplifying the variable regions 3 and 4 (V3-V4) of the 16S rRNA gene.[21] The amplicons were sequenced (300 bp paired-end) using an Illumina MiSeq sequencer at the same time in the same laboratory to avoid batch effects (Germark Biotechnology, Taichung, Taiwan). The 16S-amplicon pipeline was adapted from 16S Bacteria/Archaea SOP v1 of Microbiome Helper workflows.[22] Paired-End reAd mergeR (PEAR; version 0.9.8)[23] was used to merge paired-end reads to raw reads, then filtered low-quality reads by thresholds of sequence length ≥400 bp and quality score of 90% bases of reads ≥20. Quantitative Insight Into Microbial Ecology (QIIME; version 1.9.1) software was used to select operational taxonomic units (OTU).[24] The SILVA (version 123) 16S database[25, 26] was applied to cluster OTUs and assign taxonomy using UCLUST algorithm[27] with a 97% sequence identity threshold. Reads were dereplicated and singletons were discarded. The final OTU table was rarefied into minimum sequencing depth in the data set.
Propensity score matching
Propensity score (PS) matching[28, 29] was performed to balance confounders between the comparisons of interest (i.e., β-blocker users versus non-users) and minimize the confounding by indication resulting from nonrandom treatment study. Using a logistic regression model, β-blocker use was accessed to estimate the propensity to receive a β-blocker for each participant based on potential confounders, including age, sex, dialysis vintage, vascular access type, Bristol stool scale, dietary intake, comorbidities (diabetes mellitus, hypertension, dyslipidemia, coronary artery disease, heart failure, cerebrovascular disease, and parathyroidectomy history), concomitant drugs used (including ACEI [angiotensin converting enzyme inhibitors]/ARB [angiotensin-receptor blockers], glucose-lowering drugs [such as sulfonylurea, dipeptidyl peptidase-4 inhibitors, insulin], statin, calcium carbonate, and proton pump inhibitors), and clinical laboratory data (hemoglobin, albumin, total cholesterol, triglyceride, high sensitivity C reactive protein [hsCRP], sodium, potassium, total calcium, phosphate, parathyroid hormone, serum iron, ferritin, normalized protein catabolic rate [nPCR], and single pool Kt/V). In this study, 193 hemodialysis patients were enrolled, including 83 β-blocker users and 110 non-users (full cohort). PS-matched (1:1) analysis was used to match participants with β-blocker treatment (N=62) to participants without β-blocker treatment (N=62) (PS-matched cohort, Figure 1).
Statistical and bioinformatics analyses of microbiota
The study design is presented in Figure 1. Demographic characteristics are shown as the mean, median, or frequency, with differences between β-blocker users and non-users determined using an independent T-test or chi-squared test, as appropriate. A rarefaction curve was built to prevent methodological artifacts originating from variations in sequencing depth. The α-diversity (Shannon index and Simpson index) was estimated with the R “vegan” package and calculated using the p-value by the Kruskal-Wallis test. The β-diversity (measured by Bray-Curtis dissimilarity) was visualized through a Principal Coordinates Analysis (PCoA) and calculated using homogeneity of group dispersions by Permutational Analysis of Multivariate Dispersions (PERMDISP) to evaluate the difference/similarity of bacterial communities between the groups.[30]
Co-occurrence analysis was used to determine the relationships within communities, with core microbiome analysis performed at the genus level using MicrobiomeAnalyst[31], in which sample prevalence and relative abundance cut off values were set at 20 and 0.2%, respectively. For visualization of the internal interactions and further measurement of the microbial community, SparCC was used to calculate the Spearman correlation coefficient with the corresponding p-value between every two taxa. Microbiota community structure was assessed by co-occurrence networks built by the Sparse Correlations for Compositional data (SparCC) algorithm.[32] The p-values were estimated using a bootstrap procedure with 100 random permutations and iterations for each SparCC calculation, and correlation matrices were computed from the resampled data matrices. Only OTUs with correlation scores greater than 0.4 and p-value less than 0.05 were categorized into co-abundance groups (CAGs); these coefficients were also used to assess the length of edges on the network. An undirected network, weighted by SparCC correlation magnitude, was generated using bioinformatics tools in MicrobiomeAnalyst.[31]
The bacterial OTU difference between β-blocker users and non-users by the linear discriminant analysis (LDA) of effect size (LEfSe) analysis with more than 0.1% relative abundance and present in >30% of samples. The LEfSe analysis employed the non-parametric factorial Kruskal-Wallis test or Wilcoxon rank-sum test and LDA to identify differentially abundant taxa between the groups. Only taxa with LDA score greater than two or less than two at a p-value <0.05 were considered significantly enriched. All statistical tests are two-tailed, and a p-value < 0.05 was considered statistically significant. The random forest method[33] was performed to determine a ranked list of all bacterial taxa to identify the most predictive bacterial community to classify β-blocker users and non-users. The random forests is a supervised learning algorithm ranking OTUs based on their ability to discriminate among the groups, while accounting for the complex interrelationships in high dimensional data. The MetagenomeSeq method was also used to evaluate differential abundance in sparse marker-gene survey data using a zero-inflated Gaussian (ZIG) fit model to account for undersampling and sparsity in OTU count data after normalizing the data through cumulative sum scaling (CSS).[34] Finally, the log-transformed read counts difference of the top selected genera from the ZIG fit model between β-blocker users and non-users was analyzed in the full and PS-matched cohorts.
Co-occurrence and random forest analyses were performed by MicrobiomeAnalyst.[31] The other statistical analyses were performed using R statistical software (version 3.5.1) and STATA statistical software (version 14).