Combining amplicon sequencing and metabolomics in cirrhotic patients highlights distinctive microbiota features involved in bacterial translocation, systemic inflammation and hepatic encephalopathy

In liver cirrhosis (LC), impaired intestinal functions lead to dysbiosis and possible bacterial translocation (BT). Bacteria or their byproducts within the bloodstream can thus play a role in systemic inflammation and hepatic encephalopathy (HE). We combined 16S sequencing, NMR metabolomics and network analysis to describe the interrelationships of members of the microbiota in LC biopsies, faeces, peripheral/portal blood and faecal metabolites with clinical parameters. LC faeces and biopsies showed marked dysbiosis with a heightened proportion of Enterobacteriaceae. Our approach showed impaired faecal bacterial metabolism of short-chain fatty acids (SCFAs) and carbon/methane sources in LC, along with an enhanced stress-related response. Sixteen species, mainly belonging to the Proteobacteria phylum, were shared between LC peripheral and portal blood and were functionally linked to iron metabolism. Faecal Enterobacteriaceae and trimethylamine were positively correlated with blood proinflammatory cytokines, while Ruminococcaceae and SCFAs played a protective role. Within the peripheral blood and faeces, certain species (Stenotrophomonas pavanii, Methylobacterium extorquens) and metabolites (methanol, threonine) were positively related to HE. Cirrhotic patients thus harbour a ‘functional dysbiosis’ in the faeces and peripheral/portal blood, with specific keystone species and metabolites related to clinical markers of systemic inflammation and HE.

generated for subsequent multivariate statistical analysis. Mann-Whitney U test was used to assess significant comparisons (P ≤ 0.05).

Metagenome prediction and pathway analysis
Biom file was generated with Mothur v.1.38.1 using Greengenes database (v. 13_5_99), and used with PICRUSt 1.0.0 (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) [22] with default parameters, in order to predict Kyoto Encyclopedia of Genes and Genomes (KEGG) orthologs (KOs) from 16S V3-V4 amplicon data. Specific KOs related to significant NMR metabolites were bioinformatically assigned to each sample by means of KEGG online website (http://www.genome.jp/kegg/ko.html) and Integrated Database Retrieval System (http://www.genome.jp/dbget/), taking into account Orthology and Reactions databases for a refined search. STAMP [23] was then utilized employing twosided Welch's t-test and η 2 (effect size) to detect specific KOs having discriminant power (P ≤ 0.05) and functional relationship to NMR data of cirrhotic/controls fecal samples. Mean relative KO abundances, normalized per sample number, were computed with in-house Python scripts, and those lower than a definite threshold (2*10 -6 % for phylum, 5*10 -6 % for genus) were excluded from KO contribution graphical analysis. For peripheral and blood functional analysis, STAMP was used with Kruskal-Wallis H test, Tukey-Kramer post-hoc test (0.95), and 10% Benjamini-Hochberg FDR.

Nuclear Magnetic Resonance (NMR)
Fecal samples of controls and cirrhotics were investigated using NMR spectroscopy to solve spectra of complex mixtures and to recognize and quantify each component without chemical separation [24]. Briefly, NMR spectra of fecal samples were recorded at 27°C on a Bruker AVANCE 600 spectrometer operating at the proton frequency of 600.13 MHz and equipped with a Bruker multinuclear z-gradient inverse probehead capable of producing gradients in the z-directions with a strength of 55 G/cm. 1H spectra were referenced to methyl group signals of TSP (δ=0.00 ppm) and they were acquired by co-adding 64 transients with a recycle delay of 7 s. The residual HDO signal was suppressed using a pre-saturation. The experiment was carried out by using a 90° pulse of 11.75 μs, and 32K data points. 1H spectra were transformed with 0.5 Hz line broadening and zero filling, size 65K, manually phased, calibrated on methyl group signals of TSP (δ=0.00 ppm), and baseline corrected using the TOPSPIN v1.3 software. Spectra were prepared for statistical analysis by dividing the entire spectrum into small regions (0.02 ppm width), called "buckets". Regions with only background noise, the water resonance, and the extreme regions of spectra were excluded from the buckets. The total integral (as the sum of all 418 buckets) for each spectrum was normalized to 1000. 2D NMR experiments, namely, 1H-1H TOtal Correlation SpectroscopY (TOCSY), and 1H-13C Heteronuclear Single Quantum Coherence (HSQC) were performed using the same experimental conditions previously reported [24]. The mixing time for 1H-1H TOCSY was 80 ms. HSQC experiments were performed using a coupling constant 1JC-H of 150 Hz. Diffusion Ordered SpectroscopY (DOSY) experiment was performed using bipolar LED sequence with sine-shaped gradient of different intensities. The gradient strength was incremented in 32 steps from 2 to 95% of the maximum gradient strength (55 G/cm). The following experimental settings were used: diffusion time, Δ was 100 ms, gradient duration, δ/2 was 1.1 ms, the longitudinal eddy current delay was 25 ms, the gradient pulse recovery time was set to 100 µs. After Fourier transformation and baseline correction, the diffusion dimension was processed by means of the Bruker TOPSPIN software (version 1.3). NMR spectra were bioinformatically analyzed by in-house scripts made with Python v.2.7.11, employing probabilistic quotient normalization (PQN) [25,26], baseline removal (rollingball) and peak shifting (binning) correction. A matrix of normalized and corrected NMR peaks' areas was generated for subsequent multivariate statistical analyses

Supplementary figure S3. Hierarchical Clusterization Algorithm (HCA) of cirrhotic and control samples.
HCA was performed on all samples (x axis) taking into account the bacterial species having a mean relative abundance equal to or higher than 0.5% (y axis). Both dendrograms were generated with Bray-Curtis similarity and complete-linkage agglomeration. Colors on x axis (samples) refer to the different samples' origin, cirrhosis etiology and drug treatments (beta-blockers, diuretic, lactulose, PPI, rifaximine).

Supplementary figure S4. Beta-diversity of all samples divided by drug treatment.
PCoA analysis was performed on a distance matrix based on Bray-Curtis measure of dissimilarity, dividing fecal and peripheral blood cirrhotic samples by five drug treatments (beta-blockers, diuretic, lactulose, PPI, rifaximine). R statistic and P values were retrieved from Anosim analysis.  2C, 2D).