This study was approved by the Institutional Review Board of the First Affiliated Hospital, School of Medicine, Zhejiang University (IRB no. 2014-336). All experiments were performed in accordance with the Helsinki Declaration and the Rules of Good Clinical Practice, and no organs from executed prisoners were used in these recipients. Patients were considered for enrollment if they diagnosed with pathologically-diagnosed hepatitis B virus (HBV)-associated hepatocellular carcinoma with the Barcelona Clinic Liver Cancel (BCLC) stage A or B, and the Child -Pugh class A or B who met the Hangzhou criteria and had undergone liver transplantation more than 24 months but less than 48 months prior to the study period. Each participant filled out a baseline questionnaire (Supplementary Table S1). The inclusion criteria were: (a) aged 35–65 years with normal body weight, (b) diagnosed with hepatocellular carcinoma prior to LT and without microbial infection during the perioperative period (based on the patient’s medical history), (c) no obvious bodily discomfort, (d) no antibiotic use in the 12 weeks prior to enrollment and no probiotics and/or prebiotics after LT, and (e) FK506 tacrolimus used as the sole immunosuppressive therapy for more than 12 months. The exclusion criteria for patients were: (a) presence of severe complications such as liver abscess, recurrence of hepatocellular carcinoma (HCC), liver trauma, diabetes, hypertension, biliary and/or vascular complications through measuring serum CRP, PCT, and imaging examination, etc. (b) infection with human immunodeficiency virus, hepatitis C virus, or other types of hepatitis virus except HBV, (c) presence of any other organ-specific diseases, including intestinal diseases, pancreatic diseases, and/or tumor recurrence, (d) consumption of alcohol, tobacco, Chinese herbal medicine, and/or recreational drugs, and (e) staying up late, work fatigue and so on. Fecal samples and patient information (including data on diet, drug use, and alcohol consumption) were collected during periodic outpatient follow-up appointments between 1 December, 2014, and 31 December, 2016, at one of the two participating hospitals (First Affiliated Hospital School of Medicine, Zhejiang University, and Shulan (Hangzhou) Hospital). HCs (Group HC) matched to the age, gender, and BMI of the LT group were correspondingly screened and enrolled according to the inclusion and exclusion criteria described in our previous study (Lu et al., 2011). This study included a total of 151 fecal samples (90 from LT patients and 61 from HC) from the participants described above. Written informed consent and questionnaires addressing previous and current diseases, lifestyles, and medications were obtained from all subjects. The data of personalized perioperative medication management and professional postoperative care were extracted from electronic medical records.

Fecal samples were collected in sterile bags, refrigerated, and then taken directly to the laboratory. The samples were then divided into 200 mg aliquots, frozen rapidly in liquid nitrogen, and stored at −80°C until use.

Phenol trichloromethane was used to extract DNA from each frozen fecal sample aliquot using a bead beater to mechanically disrupt the cells, followed by phenol-chloroform extraction (Qin et al., 2014). Extracted DNA was quantified using a Qubit 2.0 Fluorometer (Invitrogen, Carlsbad, CA, United States), and molecular sizes were estimated using agarose gel electrophoresis. All fecal microbial DNA samples were diluted to a concentration of 10 ng/μl for microbial analysis.

16S rRNA gene sequences were amplified from each of the extracted DNA samples using a set of primers targeting the hyper-variable V3-V4 region (338F/806R) of the gene: 338F, 5′-barcode-ACTCCTACGGGAGGCAGCA-3′ and 806R, 5′-GGACTACHVGGGTWTCTAAT-3′. PCR amplification was performed as described in our previous study (Lu et al., 2016). DNA libraries were constructed using kits provided by Illumina Inc. according to the manufacturer’s instructions, and DNA sequencing was performed using the Illumina MiSeq 2000 platform (San Diego, CA, United States) at the State Key Laboratory for Diagnosis and Treatment of Infectious Diseases (Zhejiang University, Hangzhou, China) according to standard protocols.

Sequence data have been deposited under NCBI BioProject accession number PRJNA544155.

Clean data were extracted from raw data using USEARCH 8.0 (Edgar et al., 2011) with the following criteria: (i) sequences from each sample were extracted using each index with zero mismatch, (ii) sequences with an overlap of <50 bp were discarded, (iii) sequences in which the error rate of the overlap was >0.1 were discarded, and (iv) sequences <400 bp in length after the merge were discarded. Quality-filtered sequences were clustered into unique sequences and sorted in order of decreasing abundance to identify representative sequences using UPARSE (Edgar, 2013) according to the UPARSE operational taxonomic unit (OTU) analysis pipeline. Singletons were omitted in this step. OTUs were classified based on 97% similarity after removal of chimeric sequences using the UPARSE values (version 7.1¹ ; Zhang Y. et al., 2015). The phylogenetic affiliation of each 16S rRNA gene sequence was analyzed using RDP Classifier² (Wang et al., 2007) against the Silva (SSU132) 16S rRNA database with a confidence threshold of 70%.

Bacterial diversity was determined via sampling-based analysis of OTUs and displayed as a rarefaction curve. Bacterial richness and diversity across the samples were calculated using the following indexes: Chao1, ACE, and Shannon (Oksanen et al., 2015). To equalize the differences in sequencing depths among samples, the sequences were downsized to 5,000 per sample (20 permutations) (Zhang Q. et al., 2016). A non-parametric Mann–Whitney U-test was used to test for significant differences between two groups. Principal component analyses (PCAs) using weighted and unweighted UniFrac distance metrics were conducted, and the R package³ was used to visualize interactions among the bacterial communities from different samples (McMurdie and Holmes, 2013).

The specific characterization of fecal microbiota to distinguish taxonomic types was conducted using a linear discriminant analysis (LDA) effect size (LEfSe) method⁴ (Segata et al., 2011). Applying a normalized relative abundance matrix, LEfSe was used to: (i) identify key bacteria in fecal samples from the LT_A, LT_N, and HC groups at multiple levels in datasets, (ii) grade the key bacteria according to the results of a Mann–Whitney U-test, which determines features with significant differences in abundance between assigned taxa and uses LDA to assess the effect size of each feature (Ling et al., 2014), and (iii) visualize the results using taxonomic bar charts and cladograms. The P-values were adjusted as described by Benjamini and Hochberg (Jiang and Yu, 2017). Differences were considered significant when the false discovery rate was < 0.05.

Random forest models (Heitner et al., 2010) were introduced to identify key discriminatory OTUs between Group LT_A, LT_N, and Group HC. And to verify the key discriminatory OTUs which selected by random forest analysis, a 10-fold corss-validation analysis has been performed using rfcv function in R-package “randomForest” (R version 3.2.1). Firstly, random forest and Wilcox rank sum test were used to select differential species with both the value of Mean_decrease_in_accuracy above 0.001, and p < 0.05 by the Wilcox rank sum test (Zhang Q. et al., 2016); secondly, 10 times cross-validation analysis was performed to sift through the minimum OTU combination with the lowest error rate and the lowest number that can accurately separate the two groups; and receiver operating characteristics (ROC) analysis was then performed to measure of quality of the classification models by the R software package pROC (Robin et al., 2011). ROC curves were constructed, and the area under the curve (AUC) was used to designate the ROC effect.

Spearman’s correlation analyses were also used to assess potentially clinically relevant associations between the relative abundance of fecal bacterial genera and serum markers of liver dysfunction using Hmisc package in R.

After applying the strict inclusion and exclusion criteria described above, we enrolled 90 LT recipients (Group LT_A, n = 42; Group LT_N, n = 48) and 61 HC subjects (Group HC). There were no significant differences between the groups in terms of age, sex distribution, or body mass index. All liver recipients had previously been diagnosed as having hepatocellular carcinoma with cirrhosis, and all were HBV soluble antigen (HBsAg)-positive prior to LT. Serum levels of alanine aminotransferase, aspartate aminotransferase, and glutamyltranspeptidase were all significantly elevated in Group LT_A compared with Group HC and Group LT_N (all P < 0.001). Further details on the clinical characteristics of the subjects are provided in Table 1. Antibiotic treatments during the perioperative period (for at least 5 consecutive days) differed among patients because the use of prophylactic antimicrobial agents varied based on the perceived risks and willingness of the patients. The antimicrobial agents used included piperacillin-tazobactam, cefepime dihydrochloride, imipenem-cilastatin sodium, and micafungin sodium. The immunosuppressive therapies prescribed in the 6 months post-LT included simulect, mycophenolate mofetil, glucocorticoids, and FK506 tacrolimus; all drugs were used at doses adjusted according to patient situation based on several factors, including their patient status, platelet count, and white blood cell count. After 6 months post-LT, FK506 tacrolimus (at a serum level of 5–6 ng/ml) with or without mycophenolate mofetil (1 g/day) was used. General liver protective drugs combined with appropriately increasing dosage of anti-immunosuppressive agents were treated in those recipients of Group LT_A. And prophylaxis with high-dose hepatitis B immunoglobulin (HBIG) and Entecavir [nucleos(t)ide analogs] was used after the transplantation to suppress viral replication. As follows: for each recipient, entecavir capsule (oral medication): one capsule daily for life; and HBIG: 2000 intl units/day through intravenous drip within postoperative week 1, 2000 intl units/week through intravenous drip within postoperative week 2–5, then 800 intl units twice a week though intramuscular injection if patients fail to reach anti-HBs levels of 100 intl units/L. The whole blood concentration of HBIG after liver transplantation is needed to monitor for life. If patients fail to reach anti-HBs levels of 500 intl units/L, which tested in the 2th day after using HBIG, within the first month post-liver transplantation, dosage adjustments would be required. We enrolled recipients who were in the post-LT period of >24 months and <48 months because patients were in a stable condition during this period and were not exposed to antibiotics in the 12 weeks prior to participation in the study. In addition, all LT recipients only took FK506 tacrolimus as an immunosuppressive therapy during this period.

We generated 4,723,696 filtered sequences from the fecal samples of 90 LT recipients and 61 HCs. The qualified reads were clustered into 588 qualified species-level OTUs using 97% as the similarity cutoff. Overall, 88.3 and 60.4% of all reads could be assigned to the family and genus levels, respectively (Supplementary Dataset S1A). Species-level OTUs and species richness and diversity estimates were obtained for each microbiome (Supplementary Dataset S1B). The species richness of individual samples and of the total gut bacterial communities was estimated by rarefaction analysis. The resulting rarefaction curves indicated that the microbial richness of the sampled guts was near saturation at the applied sequencing depth (Supplementary Figures S1A,B), which was sufficient to identify most of the bacterial community members of each individual microbiome. Fecal microbiome estimated richness in ascending order was Group LT_A, Group LT_N, and Group HC (Supplementary Figure S1C). The abundance of OTUs in Group LT_A was significantly lower than that in Group LT_N and Group HC (P = 0.007 and <0.001, respectively, Figure 1A). Additionally, only 78.2% (460 bacterial OTUs) of the microbiome was shared by the three groups, with 35, 2, and 18 bacterial OTUs unique to Group HC, Group LT_A, and Group LT_N, respectively (Supplementary Figure S1D).

Subjects in Group LT_A exhibited an obviously different fecal microbiome composition compared with those in Groups LT_N and HC. In particular, the microbial diversity was significantly decreased in Group LT_A compared with the other two groups (P = 0.025 and 0.002, respectively), as estimated by the Shannon index (Figure 1B). This finding was validated by the other applied diversity parameters, namely the Chao1 and ACE indexes (Chao1 index: P < 0.000 vs. Group LT_N and P < 0.000 vs. Group HC, Figure 1C; ACE index: P = 0.002 vs. Group LT_N and P = 0.002 vs. Group HC, Figure 1D).

The Bray-Curtis (OTU number dissimilarity, Figures 2A,B), unweighted-UniFrac (qualitative, Figures 2C,D), and weighted-UniFrac (quantitative, Figures 2E,F) PCA plots, which measure the phylogenetic similarities between microbial communities, showed that the fecal microbiota of subjects in Group LT_A differed from that of the HCs, while the fecal microbiota of subjects in Group LT_N overlapped with that of subjects belonging to Groups LT_A and HC. The lower fecal microbiome diversity in Group LT compared with Group HC is also depicted in the PCA plots (Figures 2B,D,F).

Many of the taxa were differentially abundant in LT_A, LT_N, and HC. Analysis at the class level (Supplementary Figure S2A) showed that the relative abundance of Negativicutes, Gammaprotobacteria and Bacilli were significantly higher in both the LT_N and LT_A groups compared with the HC group; additionally, Fusobacteriia were significantly increased, and Firmicutes_unclassified and Lentisphaerae were significantly decreased in LT_A when compared with HC. We also performed comparisons at the phylum level (Supplementary Figure S2B) between the patient groups and healthy controls, and found that Firmicutes were significantly decreased in both the LT_N and LT_A groups compared with the HC group. Interestingly, when compared with healthy controls, LT_A showed higher relative abundance of Proteobacteria and Fusobacteria, while LT_N showed higher relative abundance of Bacteroidetes.

The results of heatmap analysis (a hierarchical clustering analysis) of the fecal microbiomes using a random forests model revealed a discriminatory intestinal microbiome between LT recipient groups and healthy subjects. Comparison of the Group LT_A and Group HC fecal microbiomes revealed 50–97%-identity OTUs assigned to 20 different families/genera that were significantly differently distributed between the two groups (Figure 3). Of these 50 OTUs, 35 OTUs corresponding to the families/genera Butyricicoccus (2 OTUs), Prevotella (2 OTUs), Clostridium_XIVb (3 OTUs), Clostridium_XIVa (1 OTU), Lachnospiraceae (12 OTUs), Faecalibacterium (1 OTU), Dorea (2 OTUs), Ruminococcaceae (5 OTUs), Romboutsia (1 OTU), Anaerostipes (1 OTU), Coprococcus (1 OTU), Blautia (3 OTUs), and Oscillibacter (1 OTU) were more abundant in Group HC than in Group LT_A. The remaining OTUs, including those corresponding to Klebsiella (1 OTU), Escherichia/Shigella (1 OTU), Clostridium_XIVa (3 OTUs), Bacteroides (4 OTUs), Fusobacterium (1 OTU), Lachnospiraceae (1 OTU), Anaerostipes (1 OTU), Erysipelotrichaceae_incertae_sedis (1 OTU), and Clostridium_XVIII (1 OTU) were more abundant in Group LT_A than in Group HC.

Comparison of the Group LT_A and Group LT_N fecal microbiomes revealed 26 most-discriminatory-OTUs. Of these, five OTUs, corresponding to Anaerostipes (1 OTU), Clostridium_IV (1 OTU), and Clostridium_XIVa (3 OTUs), were more abundant and 21 OTUs, including Prevotella (3 OTUs), Staphylococcus (1 OTU), Burkholderiales (1 OTU), Clostridium_XIVa (2 OTUs), Ruminococcaceae (2 OTUs), Lachnospiracea (8 OTUs), Bacteroidales (1 OTU), Clostridium_XVIII (1 OTU), Butyricicoccus (1 OTU), and Dorea (1 OTU), were less abundant in Group LT_A than in Group LT_N (Figure 4). Additionally, heatmap analysis of the fecal microbiomes delineated 29 distinguishing OTUs, which were assigned to 18 different genera, between groups LT_N and HC. Of these most-discriminatory-OTUs, 14 OTUs corresponding to Lachnospiraceae (7 OTUs), Dorea (1 OTU), Clostridium_XIVb (1 OTU), Ruminococcus 2 (1 OTU), Firmicutes (1 OTU), Victivallis (1 OTU), and Bacteroides (2 OTUs) were decreased, while 15 OTUs, including Megamonas (1 OTU), Prevotella (1 OTU), Lactobacillus (1 OTU), Enterococcus (1 OTU), Klebsiella (1 OTU), Veillonella (2 OTUs), Streptococccus (1 OTU), Clostridium_XIVa (1 OTU), Bacteroides (2 OTUs), Erysipelotrichaceae_incertae_sedis (1 OTU), Ruminococcaceae (1 OTU), and Lachnospiraceae (2 OTUs), were increased in the fecal microbiome of Group LT_N compared with that of Group HC (Supplementary Figure S3).

We also used LEfSe to compare the estimated phylotypes of the recipient groups with the Group HC microbiota (Figures 5A,B), with the results confirming that dysbiosis was indeed present at various phylogenetic levels. At the family level, the Group LT_A fecal microbiome was characterized by a preponderance of Bacteroidaceae, Fusobacteriaceae, Streptococcaceae, Coriobacteriaceae, and Lachnospiraceae (all LDA scores (log₁₀) >3 and P < 0.05; Supplementary Dataset S2), while the Group LT_N fecal microbiome was characterized by a higher prevalence of Enterobacteriaceae, Veillonellaceae, and Ruminococcaceae (all LDA scores (log₁₀) > 3 and P < 0.05; Supplementary Dataset S2). In comparison, the HC fecal microbiomes were dominated by Prevotellaceae, Porphyromonadaceae, Lachnospiraceae, Coriobacteriaceae, Peptostreptococcaceae, Ruminococcaceae, and Erysipelotrichaceae (all LDA scores (log₁₀) > 3 and P < 0.05; Supplementary Dataset S2). At the genus level, opportunistic pathogens (including Bacteroides, Fusobacterium, and Streptococcus) and butyrate-producing bacteria (including Clostridium_XIVa and Dorea) were enriched in the Group LT_A fecal microbiome (all LDA scores (log₁₀) > 3 and P < 0.05; Supplementary Dataset S2), while opportunistic pathogens (including Escherichia_Shigella, Klebsiella, and Veillonella), Megasphaera (lactate-utilizing bacteria), and Butyricicoccus (butyrate-producing bacteria) were enriched in the Group LT_N fecal microbiome (all LDA scores (log₁₀) > 3 and P < 0.05; Supplementary Dataset S2). In comparison, Prevotella, Collinsella, Romboutsia, Odoribacter, and butyrate-producing bacteria (including Faecalibacterium, Clostridium_IV, Ruminococcus 2, Coprococcus, Fusicatenibacter, and Clostridium_XIVb) were enriched in the HC microbiome (all LDA scores (log₁₀) > 3 and P < 0.05; Supplementary Table S2).

Additionally, Prevotella, Romboutsia, opportunistic pathogens (including Klebsiella and Escherichia_Shigella), and butyrate-producing bacteria (including Butyricicoccus, Clostridium_IV, Clostridium_XIVb, and Clostridium_XIVa) could be used to distinguish Group LT_A samples from those belonging to Group HC, with ROC-plot AUC values of 0.903 and sensitivity and specificity values of 0.863 and 0.857, respectively (Figure 6A; ROC-AUC values shown in Supplementary Dataset S2A), while Staphylococcus and Prevotella could be used to distinguish the fecal microbiomes of Group LT_A from those of Group LT_N, with ROC-AUC values of 0.801 and sensitivity and specificity values of 0.771 and 0.786, respectively (Figure 6B; ROC-AUC values shown in Supplementary Dataset S2B). Further, Victivallis, Romboutsia, Megasphaera, Megamonas, Veillonella, Lactobacillus, opportunistic pathogens (including Klebsiella, Escherichia/Shigella, Enterococcus, and Streptococcus), and butyrate-producing bacteria (including Ruminococcus 2, Fusicatenibacter, Dorea, Clostridium_XIVb, and XIVa) could be used to distinguish fecal microbiomes of Group LT_N recipients from those of Group HC, with ROC-AUC values of 0.823 and sensitivity and specificity values of 0.738 and 0.812, respectively (Figure 6C; ROC_AUC values shown in Supplementary Dataset S2C).

Next to better understand the profile of LT-associated fecal microbiome, we computed covariations between the relative abundance of the discriminatory fecal bacterial genera, and between these genera and clinical indices for liver function including ALT, AST, TB, and GGT (Supplementary Figure S4), and we noted that bacteria in most LT_A-associated bacteria including Bacteroides, Fusobacterium, Clostridium_XIVa, Streptococcus, and opportunistic pathogens (including Klebsiella and Escherichia_Shigella) especially Bacteroides and Klebsiella (Spearman’s correlation >0.3, p < 0.05 for the correlations of both genera with all serum markers; Supplementary Figure S4), positively correlated with serum markers of liver dysfunction; while HC-associated bacteria including Butyricicoccus, Clostridium_IV, Clostridium_XIVb, Prevotella, Collinsella, Romboutsia, Faecalibacterium, Ruminococcus 2, Coprococcus, Fusicatenibacter, especially Butyricicoccus and Prevotella (Spearman’s correlation <-0.3, p < 0.05 for the correlations of both genera with all serum markers; Supplementary Figure S4), negatively correlated with serum markers of liver dysfunction.