This study consecutively included 85 adult patients who underwent LTx at the Liver Transplantation Center of the First Hospital of Jilin University from 1 June 2021 to 31 March 2023. The exclusion criteria of LTx patients were as follows: 1. bloody stools; 2. severe diarrhea or constipation; 3. nontolerant to sample collection; 4. antibiotic use history within 1 month before surgery; 5. combined organ failure; 6. a history of organ transplantation other than LTx; 7. combined transplantation; 8. irregular follow-up within the first year. All LTx recipients received routine immunosuppressive therapy with tacrolimus, sirolimus, or cyclosporine. Methylprednisolone and antibiotics were routinely used after surgery. A total of 48 end-stage liver disease patients without LTx were enrolled as control (C group). All the control patients received antibiotic treatment during hospitalization. The antibiotics used in the C group were similar to those administered after LTx, especially cephalosporins. Carbapenems were used for severe infections. All the details and procedures involved in the study design were in accordance with the Declaration of Helsinki. Written informed consent was obtained from all participants before sample collection. The experimental protocol was approved by the Ethics Committee of the First Hospital of Jilin University (2021–539).

A total of 305 rectal swab samples were collected from all the subjects. Among these, 257 samples were collected from the LTx group, which were divided into four groups according to the sampling time point: before the LTx sample as the Pre group (80 samples), after LTx at discharge as the Post1 group (73 samples), the postoperative follow-up period between 1 and 3 months as the Post2 group (64 samples), and the postoperative follow-up period between 6 and 12 month as the Post3 group (40 samples). In the C group, 48 samples were collected before discharge.

Stool sample collection procedure: before collection, the perianal skin was cleaned with sterile saline and medical alcohol cotton pads, and a disposable sterile rectal swab was used to collect the sample by inserting it 5 cm into the rectum through the anus and rotating it repeatedly to ensure sufficient collection. To ensure the collection effect, each sample was collected in triplicate. Sampling was performed using standard protocols to ensure the sterile collection process. The collected samples were labeled, sealed, and immediately stored at −80°C until use.

The head of the rectal swab was repeatedly rotated and agitated in the pre-cooled saline solution. The sample was then centrifuged at ≥13,000 × g at 4°C for 10 min. The supernatant was discarded and stored on ice. DNA was extracted from the samples using the E.Z.N.A. Stool DNA Kit (Omega Bio-Tek, Inc., catalog D4015-02), according to the manufacturer’s instructions. The extracted nucleic acid concentrations were measured using the NanoDrop spectrophotometer, and the samples were stored at −80°C.

For sequencing, isolated fecal DNA was used as a template for amplification, and the V3–V4 regions of 16S rRNA were amplified by performing PCR assays using the universal bacterial primer sets 341F (5′-CCTACGGGNGGCWGCAG-3′) and 805R (5′-GACTACHVGGGTATCTAATCC-3′). The PCR reaction was performed in triplicate with a 25 μL mixture containing 25 ng of template DNA, 12.5 μL PCR PreMix, 2.5 μL of each primer, and PCR-grade water to adjust the volume. The reaction steps were as follows: 98°C for 30 s; 32 cycles of denaturation at 98°C for 10 s, annealing at 54°C for 30 s, and extension at 72°C for 45 s; and then final extension at 72°C for 10 min. The PCR products were confirmed with 2% agarose gel electrophoresis. Throughout the DNA extraction process, ultrapure water, instead of a sample solution, was used to exclude the possibility of false-positive PCR results as a negative control. The PCR products were purified using AMPure XT beads (Beckman Coulter Genomics, Danvers, MA, United States) and quantified by Qubit (Invitrogen, United States). The amplicon pools were prepared for sequencing, and the size and quantity of the amplicon library were assessed on the Agilent 2100 Bioanalyzer (Agilent, United States) and with the Library Quantification Kit for Illumina (Kapa Biosciences, Woburn, MA, United States), respectively. The libraries were sequenced on the NovaSeq 6000 PE250 platform at LC-Bio Technology Co., Ltd. (Hangzhou, China).

Analysis of paired-end data obtained by sequencing requires the following steps: 1) splitting data and removing primer and barcode sequences; 2) splicing and filtering data (cutadapt-v1.9, FLASH-v1.2.8, fqtrim-v0.94, and vsearch-v2.3.4); 3) denoising to obtain the feature sequence and abundance table (DADA2-v2019.7); 4) evaluating alpha and beta diversity (QIIME2-v2019.7); 5) species annotations (SILVA-Release138, NT-16S database, and confidence threshold for annotations: 0.7); 6) analyzing differences between comparison groups; 7) function prediction (PICRUSt-v2.2.0b, KEGG database, and MetaCyc database); 8) phenotype prediction (https://github.com/knights-lab/BugBase); 9) longitudinal analysis (q2-longitudinal); 10) correlation heatmap generation (R-3.5.2-corrplot). R (v3.5.2) packages such as ade4-v1.7.13, ggplot2-v3.2.0, and vega-v2.5.4 were mainly used for plotting the graphs. The Mann–Whitney U test was used to compare the differences between the two groups of samples with biological replicates. The Kruskal–Wallis test was used to compare the differences between the multiple groups of samples with biological replicates. Continuous correction chi-squared analysis was used to compare categorical variables. p-value < 0.05 was considered statistically significant.

In the LTx group, 83 recipients underwent DCD LTx, one underwent living donor LTx, and one underwent DCD split LTx (right half). A total of 84 recipients underwent the first transplantation, and one recipient underwent the second transplantation. The median age (range) was 50 (29–71) years, with 25 female (29.4%) and 60 male (70.6%) subjects. The median body mass index (BMI) (range) before the operation was 24.5 (16.5–36.5) kg/m². Primary diseases of all enrolled patients are as follows: 48 cases of hepatitis B, 4 cases of hepatitis C, 11 cases of alcoholic liver disease, 4 cases of drug-induced liver injury, 8 cases of primary biliary cirrhosis, 2 cases of autoimmune liver disease, and 8 cases of unknown etiology. Among the enrolled patients, 35 cases (41.2%) had liver cancer, including 29 cases of hepatitis B, 3 cases of hepatitis C, 2 cases of alcoholic liver disease, and 1 case of primary biliary cirrhosis. The patients were followed up until 31 March 2023, with an overall survival rate of 96.5%.

The control group showed comparable demographic characteristics as the LTx group: the median age (range) was 52 (27–66) years, with 14 female (29.2%) and 34 male (70.8%) subjects. The median body mass index (BMI) (range) was 23.2 (14.6–35.6) kg/m². Primary diseases were hepatitis B (23 cases), hepatitis C (1 case), alcoholic liver disease (14 cases), drug-induced liver injury (6 cases), primary biliary cirrhosis (2 cases), and unknown etiology (2 cases). Among the C group, 12 cases (25%) had liver cancer, including 9 cases of hepatitis B, 1 case of hepatitis C, and 2 cases of alcoholic liver disease.

Table 1 summarizes the baseline characteristics of these patients. Clinical information is shown in Table 2. Details of clinical measures for individual samples are provided in Supplementary Table S1.

We plotted the rarefaction curve using the number of sequences as the abscissa and the corresponding numbers of amplicon sequence variants (ASVs) as the ordinate (Supplementary Figure S1A). The rarefaction curve reflects the rationality of sequencing data and the richness of ASVs in the sample. In this study, the curves of each specimen were ended flat, indicating the sequencing was sufficient to detect rare species. The Pre group had the highest richness, Post1 group had the lowest richness, while the Post2 and Post3 groups were located in between. The rarefaction curve of the C group was higher than that of the Post1 group, while lower than that of the Post3 group (Supplementary Figure S1B).

A Venn diagram was used to show the distribution of ASVs obtained in each group (Supplementary Figure S1C). A total of 25,565 ASVs were obtained in this study, including 12,202 in the Pre group, 5,182 in the Post1 group, 7,822 in the Post2 group, 5,456 in the Post3 group, and 5,615 in the C group. In all samples, 28 phyla, 73 classes, 168 orders, 300 families, 799 genera, and 1,641 species were identified (Supplementary Table S3).

Antibiotics were commonly administered after LTx to prevent and treat postoperative infections. To distinguish the respective impacts of postoperative antibiotic usage and LTx on gut microbiota, the microbial profiles of the pre, Post1, and C groups were compared (Figure 1). Notable differences in microbial diversity were observed among the three groups (Figures 1A, B). The Shannon index was highest in the pre-group, followed by C and Post1 groups. The top 10 genera in each group are shown in Figure 1C, with the most notable change being the relative abundance of g_Enterococcus across all three groups. The top 10 taxa with statistically significant differences in relative abundance were highlighted (Figure 1D). Notably, g_Enterococcus had the highest relative abundance in the Post1 group, while g_Veillonella had the highest abundance in the C group. The other eight genera in the graph all had the highest relative abundance in the Pre group.

We compared the alpha diversity of Pre, Post1, Post2, and Post3 groups using the Chao1 index, Simpson index, and Shannon index (Figures 2A–C). Significant differences in bacterial richness and diversity were observed among the four groups (p < 0.05). The highest bacterial richness and diversity were found in the Pre group, while the lowest was found in the Post1 group, determined by all the three indexes, indicating that surgery and perioperative treatment significantly affected the microbiome structure. The alpha diversity recovered after 1–3 months subsequently but remained lower than that of the Pre group. The alpha diversity did not continue to recover in the first year to a comparable level as pre-LTx and was even significantly reduced compared to the former time point, according to the Chao1 index (Figure 2A). The Simpson index and Shannon index showed similar results as the Chao1 index, except no significant differences were found between the Post2 and Post3 groups (Figures 2B, C). Overall, the alpha diversity decreased significantly in the early postoperative period and rapidly recovered after LTx. However, the alpha diversity of the recipient’s gut microbiota cannot be fully recovered to the preoperative level in the first year after LTx.

The beta diversity differences were compared among these four groups using PCoA based on unweighted UniFrac distance (Figures 2D–H). According to PCoA among the four groups, PCo1 was 8.39% and PCo2 was 6.4% (Figure 2D). Similar to the alpha diversity, the Post1 group was segregated from other groups, suggesting the gut microbiota composition was profoundly affected. We further compared two neighboring time points’ beta diversity, and they all exhibited significant differences (Figures 2E–G), in line with the alpha diversity results. When we compared the pre-LTx gut microbiota composition with that after 6–12 months after LTx, we found that the bacterial structure was not recovered to the preoperative level (Figure 2H).

This result reflected that LTx and perioperative treatment triggered dysbiosis, which recovered shortly after transplantation; however, the diversity and composition of the gut microbiota were different from that before transplantation throughout the first year after LTx.

The cladogram shows the evolutionary tree of the top 50 microbial taxa identified in all samples (Figure 3A), where taxa that were further apart from each other on the tree diverged earlier from a common ancestor. Taxonomy abundance Circos showed the 10 dominant bacteria at the phylum level and their proportion in each group (Figure 3B). Despite the similar constitution of the top five bacteria in all groups, namely, Firmicutes, Bacteroidetes, Proteobacteria, Actinobacteria, and Fusobacteria, the ranks of Proteobacteria and Bacteroidetes were reversed in the Post1 group. Details of the proportion of each phylum in each group and the statistical differences at the level of phyla are provided in Table 3; Supplementary Figure S2.

Cluster analysis displays the similarity between groups and the relative abundance of species at the genus level in each group (Figure 3C). The left side of the figure shows the Bray–Curtis distance clustering tree structure, which then shows that the Post1 group had the lowest similarity with the other three groups, while the Pre group and the Post2 and Post3 groups shared similarity. The right side of the figure shows the relative abundance distribution of species at the genus level in each group, with the proportion of the top 10 ranked genera in each group. Among all the samples, the top 10 dominant genera were g_Enterococcus, g_Escherichia−Shigella, g_Prevotella, g_Bacteroides, g_Peptoniphilus, g_Streptococcus, g_Finegoldia, g_Faecalibacterium, g_Veillonella, and g_Klebsiella. We found that g_Enterococcus was remarkably enriched in the Post1 group. These bacteria showed a unique composition in the Post1 group.

We further demonstrated the branching evolutionary relationships of gut microbes through LEfSe analysis (Figure 3D). The linear discriminant analysis (LDA) threshold was set to 4.3. A total of 24 enriched species were identified. Among them, nine species were enriched in the Pre group, seven species were enriched in the Post1 group, five species were enriched in the Post2 group, and three species were enriched in the Post3 group.

To further elucidate the differences in species composition among groups, we conducted differential abundance testing at the genus level (Figures 4A–C). Compared with other time points, the surgery significantly triggered upregulation of g_Enterococcus, g_Actinobacteriota, and g_Ligilactobacillus, and triggered the decline of the relative abundance of most gut microbiota, including g_Prevotella, g_Peptoniphilus, g_Streptococcus, g_Finegoldia, g_Veillonella, g_Faecalibacterium, g_Anaerococcus, and g_Dialister (Figure 4A).

To investigate the difference between pre- and post-groups, we excluded the Post1 group from the analysis (Figure 4B). The top 10 differentially abundant genera among the Pre, Post2, and Post3 groups were g_Escherichia–Shigella, g_Finegoldia, g_Faecalibacterium, g_Veillonella, g_Blautia, g_Fenollaria, g_Prevotella_9, g_Bifidobacterium, g_Corynebacterium, and g_Ruminococcus_gnavus_group. The gut microbiota after 1–3 months of recovery from LTx was distinguished from counterparts of the preoperative period and 6–12 months of the postoperative period by a high proportion of g_Finegoldia, g_Veillonella, and g_Corynebacterium, and a low proportion of g_Faecalibacterium and g_Fenollaria. After 6–12 months of recovery, the feature of short recovery was disappeared as g_Finegoldia, g_Veillonella, and g_Corynebacterium decreased and g_Faecalibacterium and g_Fenollaria increased to comparable levels of the preoperative state.

When we compared the Pre group with the Post3 group (Figure 4C), we found that the top 10 differentially abundant genera were those that were highly composed in the Pre group, with an exception of g_Acidaminococcus, the proportion of which was higher in the Post3 group. We found that the bacteria that had significant changes in the Post1 group had recovered, indicating that the impact of the Post1 stage on gut microbiota was not permanent.

In this study, the ratio of statistically significant differences among the top 30 genera was 29:30, when comparing the four groups together. When the Post1 group was excluded from the analysis, this ratio became 15:30. The ratio decreased to 5:30 and 3:30 when comparing the Pre group with Post2 and Post3 groups, respectively (Table 4). With the extension of follow-up time, the gap of the gut microbiota structure between post- and pre-samples was narrowed.

In order to explore the relationship between disease severity and gut microbiota in LTx recipients, we grouped LTx patients based on their MELD-Na scores [MELD-Na = MELD +1.59 × (135 − Na)] before LTx. Patients with MELD-Na≥20 were selected in the H group, and patients with MELD-Na<20 were selected in the L group (Figure 5A). The alpha diversity was compared using the Chao1 index and Shannon index (Figures 5B, C). When comparing the high and low MELD-Na scores, there were statistical differences between H-Pre and L-Pre, but no significant differences were observed between H-Post3 and L-Post3. Moreover, when comparing the status before transplantation and at the end of the follow-up, there were statistical differences between L-Pre and L-Post3, but no statistical differences were observed between H-Pre and H-Post3. Differences in gut microbiota at the genus level between the groups were also observed (Figures 5D–G), and the most pronounced differences in gut microbiota were observed between L-Pre and H-Pre (Figure 5D). When comparing the Post3 group, except for g_Bifidobacterium, the relative abundances of differentially present species in L-Post3 and H-Post3 were generally low (Figure 5E). This suggests that the differences in gut microbiota between high and low MELD-Na patients mostly disappeared after 6–12 months. As expected, the species exhibited significant differences in both MELD groups between Pre and Post3 timepoints (Figures 5F, G), in line with the former result.

To further explore the longitudinal impact of the preoperative MELD-Na score level of LTx recipients on their postoperative gut microbiome changes, a linear mixed model was used to analyze the impact of high or low MELD-Na scores (H and L groups) on the microbial diversity and specific microorganisms. The results showed that there was no statistically significant difference in the impact of MELD-Na on the Shannon index (p = 0.472) (Figure 5H). The impact of high or low MELD-Na scores (H and L groups) on the relative abundance of top five microbial phyla and genera was analyzed. Among the genera and phyla evaluated, the MELD-Na scores only had a statistically significant effect on g_Bacteroides (p = 0.001) (Figure 5I). There was no statistical difference in other bacteria phyla and genera (Supplementary Figure S3).

The relationship between the top 10 genera of gut microbiota and concurrent clinical parameters was analyzed (Figures 6A–D). In the Pre group, g_Fenollaria was correlated with six clinical parameters: positively with CHE and Hb; negatively with TB, IDB, INR, and MELD-Na (Figure 6A). In the Post1 group, g_Enterococcus was negatively correlated with ALT, ALP, TB, and IDB (Figure 6B). The corresponding results of Post2 and Post3 groups are shown in Figures 6C, D.

In order to understand the detected bacterial phenotypes, we employed BugBase’s phenotype classification to analyze the bacterial phenotypes (Figures 7A–D). The horizontal axis in the figure represented different groups, the vertical axis represented relative abundance, and a color block area represented a phyla proportion.

As expected, the most notable difference in the bacterial phenotype was found in the Post1 group. Specifically, the relative abundance of anaerobic bacteria was significantly lower in the Post1 group than that in others (p < 0.001) (Figure 7A). The reduction of anaerobic bacteria was mainly due to the relative decrease in Firmicutes and Bacteroidetes. The anaerobic type of Actinobacteria was only present in the post1 group, which was also a feature of the gut anaerobic bacteria after LTx.

On the contrary, the relative abundance of facultatively anaerobic bacteria was significantly increased after the surgery (p < 0.001) (Figure 7B). The increase in Proteobacteria and Firmicutes contributed to this. The relative abundance of aerobic bacteria in the Post1 group was significantly promoted, which was mainly due to the increase in Firmicutes and Proteobacteria (Figure 7C). The potentially pathogenic differences of the microbiota among the groups exhibited no statistical difference among groups (Figure 7D).

In order to analyze the composition of functional genes in the existing sequenced microbial genomes, we used PICRUSt2 from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database and the MetaCyc database to perform KEGG Orthology (KO) analysis to predict the microbial metabolism function (Figures 7E, F). The figure shows the significantly different KEGG metabolic pathways at the third level and their proportions in each groups.

As expected, the Post1 group exhibited the most notable uniqueness compared to the other three groups, showing a decrease in nine pathways and an increase in six pathways (Figure 7E). The most promoted pathway was related to the fructose and mannose metabolism and phosphotransferase system pathways, suggesting that the metabolic function of gut microbiota was significantly affected after LTx.

When we analyze the metabolism difference between the Pre group and Post3 group (Figure 7F), we only found 11 metabolic pathways that were statistically different between the two groups. Then, 6–12 months after transplantation, the gut microbiota exhibited enhanced protein kinases and pathogenic Escherichia coli infection and reduced nicotinate and nicotinamide metabolism.