Patients who received a first renal transplant within 2 years and received regular surveillance for BKV were included in the study. Exclusion criteria included the following: 1) multiorgan or pediatric transplants; 2) lacking BKV detection; 3) taking high-dose antibiotics; 4) delayed graft function early after renal transplantation; 5) concurrent human immunodeficiency virus, hepatitis virus, or Mycobacterium tuberculosis infection; 6) obese patients. BKV viral load was detected by PCR in all participants. All of them received low-dose sulfamethoxazole/trimethoprim and ganciclovir for prophylaxis after transplantation. Information about posttransplant infections, rejection, and medications was collected. A total of 48 recipients were enrolled and divided into the BKV group (n = 25) and the control group (n = 23), with a ratio of proximately 1/1. Fecal samples were collected after the diagnosis of BKV infection. The study was approved by the ethics committee of Beijing Friendship Hospital (2020-P2-212-01).

Fecal samples were collected by Faeces container (Sarstedt, Germany) without preserving reagent and stored at -80°C, and processed within 2 weeks after collection.

Genome DNA was extracted from the samples using cetyltrimethylammonium bromide/sodium dodecyl sulfate (CTAB/SDS) method. DNA concentration and purity were monitored on 1% agarose gels. Extracted DNA was diluted to a concentration of 1 ng/μl using sterile water.

In this study, 16S ribosomal RNA (rRNA) genes were amplified using the specific primer Bakt_341F (5′- CCTACGGGNGGCWGCAG-3′) and Bakt_805R (5′-GACTACHVGGGTATCTAATCC-3′) (Herlemann et al., 2011). PCR reactions were performed in 30-μl reactions, containing 15 μl of Phusion^® High-Fidelity PCR Master Mix (New England Biolabs), 0.2 μl of forward and reverse primers, and 10 ng of template DNA. The thermal cycling protocol included initial denaturation at 95°C for 3 min, 25 cycles of denaturation at 95°C for 30 s, annealing at 55°C for 30 s, elongation at 72°C for 30 s, and finally 16°C for 2 min.

PCR products were mixed with an equal volume of 1× loading buffer (containing SYB green). Electrophoresis was performed on 2% agarose gel. Samples with a bright main strip around 460 bp (V3+V4) were chosen for further experiments (Vasileiadis et al., 2012).

Since we mixed the PCR products in equi-density ratios, the mixture was purified using GeneJET Gel Extraction Kit (Thermo Scientific).

NEB Next Ultra DNA Library Prep Kit for Illumina (NEB, USA) was used to generate sequencing libraries according to recommendations from the manufacturer, and index codes were added. The library quality was assessed on the Qubit@ 2.0 Fluorometer (Life Technologies, CA, USA) and Agilent Bioanalyzer 2100 system. Finally, libraries were sequenced on Illumina MiSeq platform, and 250-bp paired-end reads were generated.

The BKV infection status of each subject was confirmed by electronic medical records. Detection of BKV was performed by PCR (McNees et al., 2005). DNA was extracted from 200 µl of urine sample using the viral DNA/RNA extraction and purification kit (Xi’an Tianlong Science and Technology, Xi`an, China) according to the manufacturer’s instructions. PCR was performed using the ABI 7500 FAST Real-Time PCR System following the manufacturer’s recommendations and PCR fluorescence probe method with BKV nucleic acid detection kit (Beijing SinoMDgene Technology, Beijing, China). PCR reactions were performed in 25-μl reactions, containing 20-µl master mix and 5-µl template DNA. The thermal cycling protocol included initial denaturation at 95°C for 3 min, 40 cycles of denaturation at 94°C for 15 s, annealing, elongation, and signal acquisition at 60°C for 35 s, and finally 25°C for 1 min. BK viruria was defined as a positive result above the detectable level that the lower limit was 2.0E+03 copies/ml.

We used the random Forest (version 4.6) package of R language to build the random forest model. The sequential forward selection method was used to select the best feature set. We started with the best feature with the largest classification accurate for feature selection, and then the other feature was added one by one. Each time the model classification accuracy was evaluated, the feature with the highest accuracy was added to the model. This process was repeated until the highest accuracy for the model was achieved. The 5-fold cross-validation was used to split the data into training and test datasets in model evaluation, and the cross-validation was repeated 20 times.

Alpha diversity and beta diversity on Bray–Curtis were measured using default parameters and QIIME2 tools. Principal coordinate analysis (PCoA) was generated to reveal the divergence between groups, and analysis of similarities (ANOSIM) was applied to test the significance of the clustering based on a distance matrix of Bray–Curtis. Linear discriminant analysis (LDA) Effect Size (LEfSe) was used to identify the distinguishing bacterial taxa between groups, limiting the log LDA score >4.0.

We used SPSS (ver. 26.0, SPSS Inc., Chicago, IL, USA) for statistical analysis. Measurement data were expressed as mean [standard deviation (SD)]. Student’s t-test and chi-square test were applied to compare the quantitative variables and the categorical variables between groups, respectively. Wilcoxon rank-sum test was used to compare the abundance of microbial species between groups. Spearman correlation analysis was applied to explore the correlation between the microbiota and the CD4/CD8 ratio. A p < 0.05 was considered to be statistically significant.

The clinical characteristics of the BKV group and the control group were displayed in Table 1 . There was no significant difference in gender, age, body mass index, concurrent diabetes, medications, concurrent infections, or clinical rejection between the two groups. The BKV viral load in urine was 2.35E+10 ± 9.98E+10 copies/ml in the BKV group. The CD4/CD8 ratio was 1.18 ± 0.71 in the BKV group, lower than that in the control group (p = 0.012).

A total of 48 fecal samples were processed for 16S rRNA gene amplicon sequencing. On average, 69,683 ± 6,297 sequence reads were generated per sample in the BKV group, and 68,779 ± 6,696 sequence reads were generated per sample in the control group. There was no significant difference between the two groups (p = 0.632).

On average, the samples in the BKV group were characterized by comparable operational taxonomic unit (OTU) counts (mean: 218 ± 61 OTUs) compared with those in the control group (185 ± 70 OTUs), p = 0.083. There were a total of 390 shared OTUs between the two groups and 9 unique OTUs in the BKV group.

At the phylum level ( Figure 1A ), the abundance of Firmicutes, Bacteroidetes, Proteobacteria, Actinobacteria, and Verrucomicrobia accounted for over 99% in both groups. Firmicutes, Actinobacteria, and Verrucomicrobia were richer in the BKV group than those in the control group (63.92% vs. 58.32%, 5.08% vs. 2.69%, 3.59% vs. 2.30%), but there was no statistically significant difference (p = 0.458, p = 0.214, p = 0.469, respectively). Bacteroidetes and Proteobacteria were poorer in the BKV group than those in the control group (12.04% vs. 21.52%, 14.58% vs. 14.85%), but there was no statistically significant difference (p = 0.225 and p = 0.360, respectively). However, the Firmicutes/Bacteroidetes ratio was significantly higher in the BKV group than that in the control group (133.74 ± 306.79 vs. 26.65 ± 84.96, p = 0.046) ( Figure 1C ). The Spearman correlation analysis showed that Bacteroidetes was positively correlated with the CD4/CD8 ratio (r = 0.289, p = 0.046) ( Figure 1D ). In the top 20 dominant genera ( Figure 1B ), Romboutsia (2.26% in the BKV group vs. 0.16% in the control group, p = 0.022) and Roseburia (1.84% in the BKV group vs. 0.02% in the control group, p = 0.026) exhibited a significant difference between the two groups.

As shown in Table 2 and Figure 2 , no significant difference was observed in the microbial alpha diversity.

The PCoA plot was applied to evaluate the beta diversity based on Bray–Curtis distance analysis, which showed that the microbiota in the BKV group and the control group were clustered closely with some overlap ( Figure 3A ). The first principal component (PC1) and the second principal component (PC2) accounted for 14.64% and 12.12% of total variations, respectively. ANOSIM at the OTU level was conducted ( Figure 3B ). The R value was 0.059, while the p value was 0.028. The difference in microbiota in the inter-group was greater than that in the intra-group, indicating that the grouping was significantly meaningful.

LEfSe was used to explore the biomarkers between the BKV group and the control group. We discovered 9 distinguishing bacterial taxa between the two groups, with a log LDA score >4.0. The abundance of the class Clostridia, order Clostridiales, family Peptostreptococcaceae, Veillonellaceae, genus Romboutsia, and species uncultured bacterium of genus Romboutsia was higher, while the abundance of family Enterococcaceae, genus Enterococcus, and species uncultured bacterium of genus Enterococcus was lower in the BKV group than that in the control group. The increase and decrease of bacterial taxa abundance in the BKV group were represented by different colors ( Figure 4 ).

We established a random forest model using genus taxa to predict the BKV infectious status. The random forest model achieved the best accuracy (80.71%) to distinguish BKV-infected patients and those non-infected. The best model achieved an area under the curve (AUC) of 0.82 ( Figure 5 ). Two genera were included in the best model, which were Romboutsia and Actinomyces.

This study has limitations. Firstly, the sample size is limited. Secondly, the metagenome analysis is unused, which can provide more detailed information to explore the correlation between gut microbiota dysbiosis and BKV infection. Thirdly, a longitudinal study and metabolomics analysis are needed to explore and further understand the correlation between gut microbiota dysbiosis and the development of BKV infection.