We integrated 24 metagenomic and viral metagenomic sequencing datasets from preterm infants collected at two postnatal time points (12 samples at 14 days post-birth [Group A] and 12 samples at 28 days post-birth [Group B]; Table 1). Sequencing metadata for metagenomic and viral metagenomic datasets are summarized in Supplementary Tables S1, S2, respectively.

After removal of low-quality reads (see Methods), we obtained a total of 4.1 billion high-quality reads, with an average retention rate of 97.50% for metagenomic data and 89.52% for viral metagenomic data. Potential human DNA contamination was controlled by mapping reads to the human reference genome and discarding host-derived sequences. Following an additional round of host-filtration with BMTagger (Rotmistrovsky and Agarwala, 2011), the metagenomic datasets retained 1.7 billion host-clean reads (average 74 million reads per sample), whereas viral metagenomic datasets yielded 240 million host-clean reads (average 1 million reads per sample). Detailed statistics of host sequence removal for both datasets are provided in Supplementary Table S3.

Each sample was assembled individually to generate sample-specific contigs (see Methods). Assembly statistics for viral metagenomes, including contig number, total assembly length, N50, and maximum contig size, are summarized in Supplementary Table S4.

To characterize the taxonomic composition of the intestinal microbiome in preterm infants, we annotated metagenomic reads using Kraken and Bracken (Lu et al., 2017). In addition, to thoroughly examine the impact of viruses on premature infants, we analyzed the composition of the viral populations in both groups using CAT (von Meijenfeldt et al., 2019) based on sequencing reads (Figure 1B). The top 150 viruses, ranked by community abundance, revealed that the annotated virome was overwhelmingly composed of viruses belonging to the phylum Uroviricota, which encompasses tailed bacteriophages of the order Caudovirales. Viral relative abundances reported in this study are expressed as proportions of the total quality-controlled sequencing reads per sample, ensuring that estimates reflect the true compositional context of the gut virome.

Viral relative abundances reported in this study are expressed as proportions of the total quality-controlled sequencing reads per sample, ensuring that estimates reflect the true compositional context of the gut virome. Due to limitations in current viral reference databases, a variable fraction of viral reads could not be taxonomically assigned across samples. While exact per-sample unclassified proportions were not retained in the final analysis tables, our intermediate processing logs indicate that annotation success ranged widely—approximately 40% to over 90% of viral reads were classified in different individuals. This implies that 10% to 60% of viral sequences remained uncharacterized, likely representing novel or underrepresented viral lineages in existing databases. This implies that 10% to 60% of viral sequences remained uncharacterized, likely representing novel or underrepresented viral lineages in existing databases. Despite this variability, all identified viruses belonged to the phylum Uroviricota (100% prevalence; median relative abundance: 0.14%), which encompasses tailed bacteriophages of the order Caudovirales. All taxonomically resolved viral sequences corresponded to bacteriophages infecting bacterial hosts such as Staphylococcus, Streptococcus, Enterococcus, and Enterobacter. Key species identified included Staphylococcus phage StB20, Staphylococcus phage StB20-like, Streptococcus phage EJ-1, Enterococcus virus FL3, and Streptococcus phage YMC-2011. The phylum-level prevalence and median relative abundances reported here were computed from the taxonomic profiles in Supplementary Tables S5 (bacteria) and S6 (viruses). However, a portion of viral reads remained unidentified due to limitations in the reference databases, which primarily captured well-characterized viral species. This issue highlights the need for further expansion of viral reference databases, which will be addressed in future studies. These findings highlight the rich and varied microbial landscape within the intestinal microbiota of preterm infants, providing a foundation for future investigations into potential relationships between gut microbiota composition and clinical outcomes such as lung health.

The comparison of bacterial species composition between the 14- and 28-day groups revealed distinct distributions. The relative abundance of Staphylococcus epidermidis was significantly higher at 14 days (Paired t-test, P < 0.05), while Enterococcus faecalis was significantly higher at 28 days (Paired t-test, P < 0.05, Figure 1C; Supplementary Table S5). To provide a more transparent view of individual microbiome profiles, relative abundance profiles for each individual were displayed separately (Figures 1C, D). These individual profiles allow for a clearer understanding of the variations in microbial composition across different infants, showing that the abundance of species such as Staphylococcus epidermidis varied from negligible to over 50%, with some infants exhibiting more marked shifts in bacterial populations between the two time points. In terms of viral community relative abundance, there was a notable decrease in Staphylococcus phages from 14 to 28 days, whereas Enterococcus phages showed a general increase (Figures 1E, F; Supplementary Table S6). Streptococcus phages, which target and control Streptococcus bacteria, play a crucial role in maintaining a balanced gut microbiota by preventing bacterial overgrowth. These findings highlight the rich and varied microbial landscape within the intestinal microbiota of the study subjects, providing a robust foundation for further investigation into the role of specific microbial taxa in the health and disease states of preterm infants. This detailed microbial profiling is essential for understanding the intricate microbial ecosystems and their potential impacts on the clinical outcomes of this vulnerable population.

Analysis revealed distinct differences in bacterial and viral community composition between 14- and 28-day postnatal time points. Alpha diversity, assessed using Shannon and Simpson indices, was significantly higher in bacterial communities at 28 days compared to 14 days (Shannon index: paired t-test, P < 0.05; Figure 2A; Simpson index: paired t-test, P < 0.05; Supplementary Figure 1A). Viral communities exhibited a non-significant increasing trend in alpha diversity at the family and species levels (Shannon index: paired t-test, P > 0.05; Figure 2B; Simpson index: paired t-test, P > 0.05; Supplementary Figure 1B). Beta diversity analyses based on Bray-Curtis dissimilarity and weighted UniFrac distances did not reveal statistically significant differences between the two time points for bacterial communities (Supplementary Figure 2) or viral communities (Supplementary Figure 3). Principal Component Analysis (PCA) on species-level relative abundances highlighted structural differences in both bacterial and viral communities between groups, as evidenced by distinct clustering patterns of samples (Figures 2C, D). Overall, microbial community composition underwent significant shifts from 14 to 28 days post-birth (Supplementary Tables S6, S7).

To examine the distribution patterns and relative abundance of microbial species, we quantified shared and unique bacterial and viral species across inter-group and intra-group comparisons in this study (Figure 3). The bacterial species flower plot revealed that Group A contained 22 shared species, with approximately 41.67% of samples lacking unique species, 33.33% harboring fewer than 5 unique species, and only 25% exhibiting more than 5 unique species. In contrast, Group B exhibited 43 shared species, where 91.67% of samples contained limited unique species (Figure 3A–B; Supplementary Table S8). The analysis of viral species composition indicated that shared species outnumbered exclusive species in Group A, whereas Group B showed fewer shared species compared to unique species (Figures 3C, D; Supplementary Table S9). These findings suggest that the relative abundance of the bacterial community in the gut microbiota of preterm infants remains highly stable over time, consistent with previous observations (Maynard et al., 2012; Yatsunenko et al., 2012). The primary driver of bacterial diversity changes was identified as significant fluctuations in the relative abundance of bacterial taxa during development. Conversely, viral communities demonstrated marked variability both between and within groups, characterized by group/sample-specific viral species alongside shared members. For instance, 62.5% of viral species (145 out of 232 total) were shared between groups, while 37.5% (87 species) were group-specific. Intra-sample variations in shared and unique viral species were also evident (Figures 3E–F).

The 28-day gut-specific viral community included host-associated viruses such as Escherichia virus pro483, Megavirus chiliensis, and Spleen focus forming virus, which exhibit broad host specificity. These viruses regulate gut bacterial communities through lysis of host bacteria (e.g., Clostridium, Enterococcus, Escherichia). Additionally, bacteriophages like Clostridium phage phiCT453A and Enterobacteria phage mEp460 may facilitate horizontal transfer of resistance or metabolic genes. Giant viruses, including Megavirus chiliensis and Pandoravirus neocaledonia, potentially interact with protists or larger microbial communities, influencing gut ecosystem stability. The presence of these viruses correlates with gut microbiota remodeling and host environmental adaptability.

Among bacteriophages, six species were consistently abundant across samples, while many low-abundance species highlighted the dynamic and individualized nature of the phage community.

The 14-day gut-specific viral community predominantly targets Gram-positive bacteria (e.g., Propionibacterium, Streptococcus), suggesting roles in early microbiota establishment and stabilization. Viruses like Propionibacterium virus series may aid in initial colonization via lysis-induced nutrient release. Similarly, Pandoravirus dulcis and Pandoravirus salinus, akin to 28-day giant viruses, likely regulate gut ecosystems. The 14-day viral community exhibits high diversity, including chimeric viruses (e.g., Chimeric virus 14) and Vibrio phages (e.g., Vibrio phage ValB1MD 2), indicating complex viral dynamics during early gut development. These viruses likely support microbiota establishment and host adaptation.

Shared viral species between 14- and 28-day samples include those infecting Gram-positive (Propionibacterium, Streptococcus) and Gram-negative bacteria (Escherichia, Salmonella). Their persistent presence suggests essential roles in maintaining core gut microbiome functions. These findings align with the hypothesis that viral communities dynamically reflect host-microbe interactions, environmental exposures, or antibiotic use (Hill et al., 2017; Lugli et al., 2023).

We ranked the relative abundance of shared and unique species between groups and selected top 10 species for clustering analysis (Figure 3G–H; Supplementary Figures 6–7; Supplementary Tables S9–S10). Microbial abundance differences between groups indicated relative stability in bacterial and viral communities at these time points, though denser sampling is required to confirm longitudinal trends. However, notable shifts in abundance structures likely contributed to diversity variations between groups. These results corroborate the view that early-life gut microbes gradually respond to internal/external influences, altering bacterial/viral community prevalence (de Muinck and Trosvik, 2018).

The gut microbiota of preterm infants undergo dynamic changes during early development, consistent with previous studies documenting rapid microbial succession in neonatal gut ecosystems (Henderickx et al., 2019; Toubon et al., 2022; Westaway et al., 2022). To investigate these changes, we analyzed the relative abundance and diversity of bacterial and viral species in the gut microbiota of preterm infants at 14- and 28-day postnatal time points (Figure 4). Notable shifts were observed, including a reduction in Staphylococcus species prevalence, particularly a marked decline in Staphylococcus epidermidis relative abundance (Figure 4A; Supplementary Table S12), and a significant increase in Enterococcus faecalis abundance at 28 days compared to 14 days. These shifts may be attributed to factors such as antibiotic therapy (Table 1) and potential probiotic supplementation. For instance, Lactobacillus rhamnosus-based probiotics have been shown to enhance microbial diversity while reducing Enterobacteriaceae and Staphylococcaceae abundance, thereby promoting a healthier gut microbiota (Martí et al., 2021). The observed increase in Enterococcus faecalis may reflect its role as an opportunistic colonizer in the preterm gut. Given the underdeveloped microbiota and immune systems, its presence could represent a natural phase of microbial succession. However, further investigation into its functional impact on gut health and immune development is warranted.

Regarding viral communities, we observed changes in the relative abundance of Staphylococcus- and Enterococcus-related bacteriophages between 14 and 28 days (Figure 4B; Supplementary Table S13), reflecting fluctuations in their bacterial host populations. These findings underscore the dynamic phage-bacteria interactions within the gut microbiota. Although our analysis did not explicitly perform correlation tests between time and microbial shifts, paired t-tests confirmed statistically significant (P < 0.05) temporal changes in viral and bacterial abundance. These results highlight the age-dependent evolution of microbial and viral communities in preterm infants. We included a broad range of bacterial species, including low-abundance Staphylococcus species (excluding S. epidermidis) (Figure 4), to comprehensively capture microbial diversity and colonization dynamics. While these species may contribute minimally to overall microbiota composition, their inclusion provides a more complete picture of the early gut ecosystem, potentially revealing minor contributors critical to initial colonization. Additionally, even low-abundance species may play ecological roles or influence microbial balance, particularly in the context of immature immune systems in preterm infants. The proportion of unidentified viral reads in our dataset is detailed in Supplementary Table S7. Although most viral species were classified, a fraction of reads remained unassigned due to limitations in reference databases, emphasizing the need for expanded viral reference datasets in future studies to improve identification accuracy.

To further explore functional differences between bacterial and viral communities, we conducted KEGG pathway enrichment analysis (Figure 5). Among 21 significantly enriched pathways, those related to the phosphotransferase system (PTS) and starch/sucrose metabolism were prominent at 28 days, indicating adaptation to dietary carbohydrate utilization. In contrast, the Staphylococcus aureus infection pathway remained consistently low across both time points (Figures 5A, B; Supplementary Tables S13–S14). Statistical significance was determined using FDR (False Discovery Rate) correction, which controls false positives during multiple comparisons, ensuring reliable identification of significant pathways. Functional analysis of carbohydrate-active enzymes (CAZymes) revealed age-dependent shifts in microbial metabolic activity: Auxiliary Activities (AA), GlycosylTransferases (GT), Carbohydrate Esterases (CE), and Carbohydrate-binding Modules (CBM) predominated at 14 days, while Glycoside Hydrolases (GH) and Polysaccharide Lyases (PL) became more prevalent by 28 days (Figure 5C; Supplementary Table S15). These findings suggest developmental changes in carbohydrate metabolism, potentially influencing complex carbohydrate degradation and microbial ecology in the maturing gut (Loss et al., 2023).

To better characterize functional differences in viral community relative abundance between 14- and 28-day postnatal time points in preterm infants, we applied Linear Discriminant Analysis Effect Size (LEfSe) to identify significant functional biomarkers associated with microbial shifts. These biomarkers may offer insights into viral influences on gut health, though further research is needed to clarify their clinical relevance in preterm populations. In our study, LEfSe was used to analyze viral functional relative abundance differences in preterm infant gut microbiota across these two time points. The analysis began with non-parametric tests (e.g., Kruskal-Wallis rank-sum test) on viral functional data, followed by Linear Discriminant Analysis (LDA) to quantify effect sizes and compute LDA scores (Figure 6A; Supplementary Table S16). LEfSe results revealed that antimicrobial resistance and cationic antimicrobial peptide (CAMP) resistance pathways were significantly enriched at 14 days, while starch/sucrose metabolism, PTS, carbohydrate metabolism, and energy metabolism pathways dominated at 28 days. These findings correlated with bacterial community functional trends, suggesting auxiliary metabolism genes (AMGs) may mediate these shifts, consistent with bacterial functional relative abundance differences observed between groups. The dynamics of bacteriophage infectivity are critical for understanding microbial ecosystem stability. We examined phage infectivity patterns at 14- and 28-day postnatal periods and found a higher proportion of temperate phages at 28 days compared to 14 days. This suggests a temporal increase in temperate phage prevalence during this window. However, due to sampling limited to two time points, it remains unclear whether temperate phage proportions would continue to rise or decline beyond 28 days. Longitudinal studies with additional time points (e.g., 36 days or later) are necessary to establish long-term phage infectivity trends (Figure 6B). While low viral loads and limited sample sizes precluded statistically significant p-values, median value differences between groups were evident.

Bacteriophages, as viruses that replicate within bacterial hosts, often influence host metabolism through auxiliary metabolism genes (AMGs). These AMGs modulate bacterial metabolic pathways to support phage replication and propagation. We classified AMGs from both time points based on their involvement in metabolic processes (Figures 6C, D). Both groups showed dominant AMG activity in amino acid metabolism, including genes encoding enzymes for amino acid synthesis, degradation, or transport. At 28 days, carbohydrate, energy metabolism, and cofactor/vitamin pathways exhibited significant increases. While phage-related quorum sensing or cancer genes were not explicitly identified, our results suggest phages may alter host bacterial metabolism to facilitate replication and reshape the microbial ecosystem.

Notably, nucleotide metabolism and aromatic compound metabolism emerged in 28-day postnatal phage AMGs. Nucleotide metabolism genes (e.g., synthesis, degradation, or salvage pathways) and aromatic compound-metabolizing AMGs encode enzymes enabling phages to utilize, degrade, or modify aromatic substrates during host infection. In contrast, terpenoid/polyketide metabolism was prominent in 14-day postnatal phage AMGs. Although some AMGs encode terpenoid-degrading enzymes, terpenoids primarily function as secondary metabolites supporting microbial survival (e.g., antimicrobial defense, signaling) rather than directly fueling phage replication. The observed terpenoid decrease at 28 days likely reflects indirect effects on phage activity, such as ecosystem-level shifts or bacterial host availability, aligning with increased temperate phage prevalence at 28 days. Previous studies have shown that terpenoids regulate microbial interactions through antimicrobial and ecological effects (Dorman and Deans, 2000; Gershenzon and Dudareva, 2007; Zhou et al., 2024). Additionally, terpenoid-degrading AMGs may modulate host microbiota rather than directly provide phage energy (Hurwitz and U’Ren, 2016). More strikingly, most phage AMG categories showed increased abundance at 28 days compared to 14 days. These findings highlight dynamic changes in viral functional profiles, phage infectivity, and AMG activity in preterm infant gut microbiota during the first month of life. The observed shifts reflect adaptive responses to evolving microbial environments and metabolic demands of both viral and bacterial communities.

Predicting interactions between phages and clinical variables is crucial for advancing our understanding of microbial dynamics and their influence on health outcomes. Phages modulate microbial communities by specifically infecting and lysing bacterial hosts. To explore these interactions, we constructed a microbial interaction network linking clinical data from preterm infants with bacterial species, including 782 shared species across groups (Figure 7A; Supplementary Table S17). The Enterobacter roggenkampii network exhibited the strongest associations with daily milk intake (14–28 days post-birth) and oral feeding initiation timing. However, only 4 of the 12 infants had documented initiation times within this period, limiting the statistical robustness of this relationship. Notably, 7 out of 12 infants were twins, and microbiome data from twins may exhibit non-independence due to shared genetic and environmental factors, potentially introducing bias into statistical analyses. Given the small sample size and use of traditional statistical methods that neglect twin-related correlations, the impact of twin pairing remains a minor limitation. Future studies with larger cohorts, more comprehensive clinical datasets, and statistical frameworks accounting for non-independent twin data (e.g., mixed-effects models) are needed to validate these associations and clarify the role of oral feeding in shaping microbial dynamics.

As anticipated, antibiotic use within four weeks post-birth showed strong connections with Pseudomonas aeruginosa, a well-documented opportunistic pathogen linked to hospital-acquired infections and antibiotic resistance. Network analysis revealed significant associations between microbial species and clinical parameters: Streptococcus thermophilus demonstrated the strongest connectivity with post-four-week antibiotic duration and types, while S. pneumoniae exhibited weaker yet statistically significant correlations with admission temperature (Figure 7B). Furthermore, the heatmap of top 20 species highlighted key bacterial taxa associated with clinical variables, including milk intake and time to regain birth weight (Figure 7C). This analysis illustrates the complex interplay between microbial communities and clinical outcomes in preterm infants. However, S. pneumoniae was not among the top 20 species ranked by the absolute sum of correlation coefficients across all clinical variables (Figure 7C). These findings underscore the importance of elucidating phage-bacteria interactions and their clinical relevance to inform therapeutic strategies for preterm infants.

The collected fecal samples were immediately sent to the laboratory on the same floor. The stool was divided into 2–4 tubes based on the amount collected, placed in sterilized Eppendorf tubes, and then stored in the laboratory –80°C refrigerator for further Uniform DNA extraction of fecal microbiota. This study was approved by the Shenzhen People’s Hospital Institutional Review Board (LL-KY 2022288).

For metagenome library construction, 1 µg of extracted genomic DNA was fragmented to an average size of ~250 bp using a Bioruptor sonicator (Diagenode). The fragmented DNA was then subjected to end repair, blunting, and phosphorylation using T4 DNA polymerase, Klenow Fragment, and T4 polynucleotide kinase. The DNA fragments were then 3’ adenylated using the Klenow Fragment (3’-5’ exo-) and subsequently ligated to adaptors using T4 DNA Ligase. After each reaction step, the DNA was purified using the QIAquick PCR purification kit (Qiagen).The final libraries were generated by PCR amplification using adapter-compatible barcode primers. For Metavirus, samples were resuspended in 10ml of SM buffer and homogenized. After 5min chilling on ice, the samples were centrifuged at 5000rpm for 10min at 4°C. The supernatants were transferred to a new tube and centrifuged again. The supernatants were then filtered twice with a 0.45um filters. A final concentration of 0.5M NaCl and 10% of PEG8000 were then added, and incubated overnight at 4°C. Then the samples were centrifuged at 5000rpm for 20min at 4°C, collecting the precipitate. Pellets were resuspended in 400ul of SM buffer and the same volume of chloroform. After emulsification, the samples were centrifuged at 2500g for 5min. The aqueous phase was transferred into a new tube and mixed with DNase and RNase incubating at 37°C for 1hour. VLP nucleic acid was then extracted using QIAamp Viral RNA Mini Kit according to the manufacture’ s instructions. After reverse transcription for RNA nucleic acid, both DNA and cDNA were amplified using MDA technology, respectively. The amplified dsDNA was then used for library construction by E-GENE Tech Co., Ltd. Both libraries were quantified with an Agilent 2100 Bioanalyzer (Agilent Technologies) and real-time PCR assays. Sequencing was performed at E-GENE in Shenzhen using the Illumina sequencing platform.

Based on unique barcode and primer sequences, paired-end reads were accurately mapped to their respective samples. High-quality clean reads were generated by removing adapter sequences and low-quality reads using Trimmomatic (Bolger et al., 2014, v0.38)(V0.38). with the following parameters: SLIDINGWINDOW:5:15, HEADCROP:3, AVGQUAL:15, LEADING:5, TRAILING:5, MINLEN:80. Sequencing was performed on the Illumina NovaSeq 6000 platform using 150-bp paired-end reads, employing Illumina’s latest S4 flow cell chemistry and reagent kits to ensure high throughput and data accuracy. To monitor contamination, kit-negative controls (no-template controls) and sequencing-negative controls (blank library preparations) were processed in parallel with experimental samples during DNA extraction, library construction, and sequencing. Analysis of control data revealed negligible read counts, confirming minimal contamination risk. For metagenomic sequencing data, clean reads were assembled into contigs using MEGAHIT (v1.0.6) (Li et al., 2016).Host DNA was removed via the NCBI Human Sequence Removal protocol using a human reference genome (hg38). Taxonomic annotation was performed using Kraken2 & Bracken with the Kraken2 Virus database. Metavirome data were processed using identical methods: host DNA was filtered using the NCBI Human Sequence Removal protocol with hg38 reference, reads were assembled into contigs via MEGAHIT (v1.0.6), and taxonomic classification was conducted using CAT (v5.3.2) through the LCA (lowest common ancestor) algorithm.

Alpha and beta diversity analyses were conducted to assess microbial shifts between 14 and 28 days post-birth. Alpha diversity, reflecting species richness and evenness within each sample, was calculated using the Shannon index. Statistical comparisons of alpha diversity between the two time points were performed using paired t-tests. Beta diversity was analyzed using two complementary approaches: Bray-Curtis distance, which measures compositional dissimilarity based on relative abundance, and weighted UniFrac distance, which accounts for phylogenetic relationships among taxa. Differences in beta diversity were assessed using ANOSIM (Bray-Curtis distance) and PERMANOVA (Permutational Multivariate Analysis of Variance). ANOSIM provided an R-value (ranging from −1 to +1),where 0 suggests the null hypothesis cannot be rejected (Clarke, 1993)), where 0 suggests no difference between groups, and PERMANOVA provided an R² value (ranging from 0 to 1), quantifying the proportion of variation explained by group differences. Sequencing depth normalization was conducted using rarefaction, and data were transformed using centered log-ratio (CLR) transformation for beta diversity analysis. Differences in community structure based on beta diversity were visualized using principal coordinate analysis (PCoA) with the ape package and non-metric multidimensional scaling (NMDS) with the VEGAN package(v2.5-3) (Oksanen et al., 2015), both using Bray-Curtis distance. Statistical significance for both PCoA and NMDS was assessed using PERMANOVA. Phylum, genus, and species differences between groups were identified using STAMP (v2.1.3) (Parks et al., 2014), with Welch’s t-test used for statistical comparisons. All analyses defined statistical significance as P < 0.05.

Open reading frames (ORFs) were predicted by prodigal (v2.6.3) form contigs and the clean reads were aligned to ORFs to calculate the mapped read count. Thereafter, the ORFs were clustered to remove redundancy for building a set of representative genes through CD-HIT software (Fu et al., 2012) (v4.7). To avoid bias related to variation in the ORF size, both ORF sequence length and sequencing depth were included in the data normalization process prior to statistical comparisons (Hu et al., 2020). The normalized abundance matrix (G) was calculated for all samples:

where m is the number of samples, n is the number of representative genes, gkh is the normalized abundance of representative gene h for sample k, p is the total number of ORFs clustered for representative gene h and sample k, Dki is the mapped reads count of ORF i for sample k, and Lki is the length of ORF i for sample k.

ORFs were annotated with functional information using the Annotation softwarer KofamKOALA of Kyoto Encyclopedia of Genes and Genomes (KEGG) (Aramaki et al., 2020) based on the KEGG database with default parameters. COG annotation of ORFs was performed using eggnog-emapper software. ARGs were characterized by mapping ORFs to the Comprehensive Antibiotic Resistance Database (CARD) (V3.0.0) using the Resistance Gene Identifier (RGI) application. The CARD is a rigorously curated collection of characterized, peer-reviewed antibiotic resistance determinants (McArthur et al., 2013). Diamond was used to map the ORFs to the MGE database. The MGE database contains genes with 278 different gene name annotations (annotated as IS*, ISCR*, intI1, int2, istA*, istB*, qacEdelta, tniA*, tniB*, tnpA*, and Tn916 transposon in the NCBI nucleotide database) and more than 2000 unique sequences, excluding the sequences from the PlasmidFinder database (Pärnänen et al., 2018). Subsets of ORFs that were mapped to MGEs were used to produce the MGE profiles. The analysis of Metagenome also refer to this method.

Spearman correlation coefficients were calculated to assess the relationships between microbial species and clinical variables. Correlation networks were constructed focusing on interactions with |R| > 0.7 and P < 0.05.

The correlation analysis between virus and clinical data was performed using R package corrplot.