To understand the role of the urinary virome in renal transplant and BK disease, we performed metagenomic sequencing to profile DNA viral communities of 65 urine specimens from 23 adult kidney transplant recipients (Missouri, USA). This study is a retrospective analysis of urine sample collected between 2000–2002 from kidney transplant recipient patients enrolled in a prospective study to prevent BK nephropathy (16). Eleven individuals were clinically diagnosed with BK viremia (BKV+), the hallmark of progression towards BK virus-associated nephropathy (BKVAN), and 12 recipients without BK viremia (BKV-). Fifteen transplant recipients were male (72.7% of BK viremia recipients) and 8 were female (27.3% of BK viremia recipients). The median age in this study was 57 years, 3 African American patients (13%) and 20 White patients (87.0%) (Fig. 1A). Repeated urine specimens were collected from transplant recipients between 0–14 weeks post-transplantation to assess changes over time (Fig. 1B). The most abundant viral contigs were from viral families Polyomoviridae (31%), Anelloviridae (10%), Circoviridae (9%) and Peduoviridae (4). There were also a large proportion of the contigs that were unclassified bacteriophages from class Caudoviricetes (28%) or were unclassified viruses (9%).
Urinary virome signatures associated with BK Viremia
We first wanted to assess if there were interactions between the metadata of patient sex with BKV status, patient sex with post-transplant time, and BKV status with time for virome richness, alpha and beta diversity. We applied LME (richness and alpha diversity) and PERMANOVA (beta diversity) models but did not find statistically significant interactions between patient sex and BKV status or time for virome richness, alpha or beta diversity (p-values of patient sex*BKV/patient sex*time; richness p = 0.24/0.57, alpha diversity p = 0.96/0.065, beta diversity p = 0.096/0.60). We found significant interactions between BKV status and time for richness and alpha diversity, but not beta diversity (richness p = 0.012, alpha diversity p < 2.28E-05, beta diversity p = 0.50; Fig. 2A-B). Since the models indicate that there is an interaction between BKV status and post-transplant time, we stratified data by BKV status (BKV + and BKV-) for richness and alpha diversity. When stratified, we found the richness in the BKV + to be significantly increasing over time, however the BKV- richness did not change throughout the 14 weeks (BKV + p = 0.0059, BKV- p = 0.932; Fig. 2A). Instead, the richness in the BKV- samples was significantly different by patient sex, with the females having more richness in urine than males (p = 0.0063; Fig. 2C). Interestingly, virome alpha diversity increased significantly over time in BKV + patients, but this was reversed in BKV- patients (BKV + p = 0.0009, BKV- p = 0.014; Fig. 2B). Virome beta diversity clustered distinctly by BKV status (p = 0.001; Fig. 2D). Taken together, this suggests that the substantial changes in the urinary virome are associated with BKV in renal transplant recipients.
Since BK polyomavirus infection itself is linked to the disease, we sought to control for this by removing BK Polyomavirus (i.e., BKV dominated signals) and rerun the LME and PERMANOVA models above. By excluding BKV, this would ascertain whether there were other bona fide urinary virome signatures. Beta diversity models showed that there was an interaction between patient sex and BKV status (p = 0.033). When stratified by patient sex, we found that the virome beta diversity differed significantly by BKV status in males but not females (female BKV status p = 0.834, male BKV status p = 0.008; Fig. 3A). Similarly, virome richness differed by BKV status as well as by patient sex. The urinary virome richness in females was higher than males, and was dependent on BKV status (patient sex p = 0.0023, BKV status p = 0.023; Fig. 3B-C). Virome alpha diversity was still significant for the interaction between BKV status and time. The opposing trajectories in virome alpha diversity was still consistent, showing that BKV- significantly decreasing over time (p = 0.013) and BKV + increasing over time though this was no longer statistically significant (p = 0.31) (Fig. 3D). Overall, we found significant differences in the virome by BKV status both with and without data dominated in BK polyomavirus. This confirms that the BKV is a strong influence on the virome of the kidney transplant patients.
Urinary virome community status
We next characterized the urinary virome by community state profiles. Given the BKV dominance in some samples, we clustered them separately from the rest. We found that BK dominant viromes clustered into three distinct clusters that were made up BK polyomavirus contigs specific to their cluster (Fig. 4A). Since there are four main subtypes of BK polyomavirus (subtypes I, II, III and IV), we hypothesized that the clusters were driven by BKV subtypes. To test this, we performed a phylogenetic analysis on the VP1 region which can be used to classify subtypes from our samples and reference sequences. We found that community state group 1 samples were comprised of BK subtype IB2, group 2 contigs were from subtype IB1 and group three contigs were from subtype IA. Further, the BKV sequences from same individual at multiple different times were phylogenetically more closely-related to one another consistent with a chronic BK infection (Fig. 4B).
In the samples without BK dominant samples, we found 6 distinct urinary community states (Fig. 4C, Supplementary Fig. 1A). Group 1 was most abundant in 5 Caudoviricetes contigs (between 13 − 3%, sum 29%), 2 Anelloviridae contigs (10% and 4%, sum 14%) 2 Autographiviridae contigs (both 8%, sum 16%) and 1 Peduoviridae contig (5%). Group 2 was highly abundant and prevalent in 10 Circoviridae contigs (5% − 3%, sum 33%), one unclassified contig (5%), one Inoviridae contig (4%), and one Rountreeviridae contig (3%). Group 3 was dominated by two Anelloviridae contigs (48%, and 42%, sum 91%). Group 4 was abundant in one Peduoviridae (5%), one Microviridae contigs (4%), two Caudoviricetes contigs (both 3%, sum 6%), three Anelloviridae contigs (all 3%, sum 9%), and one Anaerodiviridae contig (3%). Group 5 was abundant in two Anelloviridae contigs (35% and 17%, sum 53%), different from group 3 Anelloviridae contigs. Group 6 was dominate one unclassified contig (92%). We found no significant association between groups and BKV status, patient sex, or time by multinomial logit models.
Plasma anellovirus loads have been implicated in immunosuppression and allograft rejection in solid organ transplant recipients (17–19). We found that the counts of anelloviridae reads in urine was significantly higher in BKV- specimens compared to BKV + in the sequencing dataset (p = 0.0017) (Supplementary Fig. 1B). We also found the females to have significantly more anelloviridae in their urine in comparison to males (p = 0.0038, Supplementary Fig. 1C). Using qPCR assays to quantify anellovirus load in the urine specimens, we detected urine anellovirus load at higher median levels in BKV- specimens (median 15,786 copies/ml) than compared to BKV + specimens (median 981.9 copies/ml). However, the difference was not statistically significant (p = 0.149) (Supplementary Fig. 1D). When comparing the anellovirus qPCR load by patient sex, we found the females (median 28,523 copies/ml) to have significantly higher anellovirus than males (117.8 copies/ml, p = 0.0002, Supplementary Fig. 1E).
Because anellovirus contigs were found in many samples with no individual cluster associated with all contigs, we decided to create a phylogenetic tree of the anelloviruses from our samples. We found that anellovirus species TTV5, TTV13, TTV16, TTV28, and TTV29 were most closely related to the anellovirus contigs from our samples. We also found that two patients (18, 13), both females from control group, had multiple analogisms species throughout their time points in this study (Fig. 4D).
Differential signatures of the urinary virome
We next wanted to determine if there were specific viruses associated with BKV status, patient sex, or time. When looking all the data together, we found 4 Polyomavirus contigs to be the most important factor differentiating by BKV status (p-values between 0.0002–0.001, Fig. 5A). We also found one unclassified contig to be important by time, with it being high in abundance and prevalence in samples before transplant and then decreasing over the next few weeks following transplant and then seeming to recover at week 12 onwards (p = 0.003, Fig. 5B). We also tried to run masslin2 on data stratified by BKV dominant and non- dominant, but we did not find anything significant in either group by BKV status, patient sex, or time. We then stratified the data by patient sex and found no differential contigs by BKV status, patient sex, or time in males. However, 11 Polyomaviridae contigs were found to be differentiating the females by BKV status (p-values between 0.0003–0.009, Fig. 5C). We also found 7 Circoviridae contigs (p-values between 0.0031–0.004) and 1 Anelloviridae contig (p = 0.02) differentiating the females by time (Fig. 5D). Circoviridae contigs were high in abundance and prevalence pre-transplant and then started to decrease over the weeks following transplant and then finally seem to be recovering towards week 14. Taken together, Polyomaviridae contigs differentiated samples by BKV stats as expected, and we also found contigs, mainly Circoviridae, which seem to be decrease following transplant and recover towards week 14.