This study includes five metagenomic datasets, including IBS datasets sequenced from recruited subjects and four other datasets from the public database. The dataset of 22 IBS patients and 15 healthy controls have been described in our previous paper (Xu et al., 2020). Briefly, these subjects were recruited at the Outpatient Department of Gastroenterology of Peking University Third Hospital. The studies involving human participants were reviewed and approved by the Ethics Committee of Peking University People’s Hospital (No.2017PHB105-01). The participants provided their written informed consent to participate in this study. IBS patients should meet the standard of the Rome III criteria. Exclusion criteria included organic gastrointestinal or systemic diseases, use of antibiotics or antidepressants within a month, and use of probiotics, laxatives, or antidiarrheal drugs for more than 3 days during the previous 2 weeks. Among these individuals, nine cases and nine controls have corresponding metabolomic data (non-targeted metabolomics profiling on serums). Details of the sample collection, DNA sequencing, and metabolomic assay are documented in Supplementary Methods in Supplementary Materials.

Another four metagenomic published datasets were also obtained from studies related to the following diseases: T2D (Qin et al., 2012), CD (Lewis et al., 2015), CRC, and LC (Qin et al., 2014; Yu et al., 2017). The selection of these four diseases was mainly in consideration of metabolic or bowel dysfunctional diseases and data availability when we performed the research. Sample information is included in Table S4. Except for one CD dataset which is an American cohort, all are Chinese cohorts.

We began with the raw reads and processed them in a uniform pipeline. Reads were first quality controlled by prinseq-lite, with arguments -ns_max_p 10 and -min_qual_mean 25. Then, human sequences were removed by mapping reads to human reference genome GRCh38 with bowtie2 using the argument -very-fast (Langmead and Salzberg, 2012). The remaining reads were considered as clean reads. To assure the quality of these samples, only samples of clean reads fastq files larger than 2 Gb and contig N50 lengths longer than 1 kb were included, which resulted in 27 cases versus 31 controls in the T2D dataset, 148 cases versus 18 controls in the CD dataset, 52 cases versus 51 controls in the CRC dataset, and 40 cases versus 34 controls in the LC dataset. The 37 sequenced samples in the IBS dataset all met the sample quality standard.

The overall analysis pipeline is shown in Figure 1A. After quality control and removing host sequences that mapped to the human genome, clean reads’ file size ranged from 2 to 28 GB. The metagenomic reads were then assembled into contigs by metaSPAdes (contigs longer than 1.5 kb were kept), and genes were predicted from these contigs by MetaGeneMark (Zhu et al., 2010; Nurk et al., 2017). We assumed that gene abundance was approximately proportional to the actual organism abundance. Thus, these genes were used to estimate the relative abundance of taxon among viruses and bacteria, respectively. For viral taxonomic annotation, genes were first aligned against the NCBI RefSeq non-redundant protein database (O'Leary et al., 2016) and then to the hidden Markov models of viral protein families generated from the JGI Earth’s virome project (Paez-Espino et al., 2016), in which 167,042 protein-coding genes from 2,353 isolated viral genomes were clustered into 14,296 protein families. To remove false positives that originated from bacteria, we removed genes that had significant hits to human gut bacteria contigs published in Forster’s study (Forster et al., 2019) (blastn, with arguments: -perc_identity 98, -qcov_hsp_perc 90, -evalue 1e-10, -max_target_seqs 1). The remaining union of genes that had hits against reference viral proteins and viral protein families were considered viral genes. We benchmarked the viral identification process by a simulated dataset, which revealed a specificity (true positive rate in the predicted viral genes) of 98.8% and a recall rate (true positive rate in all viral genes) of 54.0% (see details of benchmark in Supplementary Methods). Meanwhile, the non-viral genes were annotated with Kaiju (Menzel et al., 2016), which maps sequences to the NCBI RefSeq protein database containing bacteria and archaea proteins.

All our taxonomic classifications were based on predicted genes. We first evaluated the contig abundances by mapping reads to contigs and calculated reads per kilobase: a i = x i L i × 1000 , where a_i is the abundance of contig i in sample S, x_i is the number of reads mapped to contig i, and L_i is the length of contig i. The gene abundance was calculated after adjusting the corresponding contig abundance by gene number in that contig: g i j = a i n i , where n_i is the number of genes in contig i. The abundance of a taxonomic unit was then added by the abundance of genes with the same taxonomic annotations: t = Σ g_{i j} , t ' = t ∑ t , where t′ is the relative abundance of a taxonomy unit in sample S.

Based on the abundance matrix constructed by calculating the abundance of each feature, we used SparCC to measure the correlation between each pair of features (Friedman and Alm, 2012). One of the major factors that we need to take into consideration when inferring co-occurrence relationships of the microbiome is composition bias, which refers to the bias caused by the relative abundances that sum to 1 so that fractions tend to be negatively correlated regardless of their true relationships. SparCC is designed to estimate correlations of compositional data by the log transformation of pairwise data. Since SparCC adds a small value to each zero value to perform log transformation, correlations for elements that appear only in a few samples may not be reliable, so we removed elements that appear in less than 20% of the sample within a group.

After constructing networks of each case and control group by SparCC, the differential interaction between case and control groups of a cohort was divided into two parts: links that existed only in the network of the case group and links that existed in both networks but had significantly different correlation coefficients [r(s_j ,s_j ′)]. To measure the significance of the difference of r, we calculated the difference (D, absolute value) of r between the same links that existed both in the network of case and control and compared it with the bootstrapped difference of r between cases and 100 shuffled co-abundance matrixes generated in SparCC. If D was beyond the 95% quantile of the bootstrapped difference, the interaction represented by this link was significantly different between cases and controls. The case-specific networks selected in this way were further filtered with SpiecEasi to get more reliable links and reduce indirect interactions (Kurtz et al., 2015).

We inferred the pan- and core-networks at the family level, since the pool of the genus-level networks had too many links to be visualized. The pan- networks were constructed by the pool of the links of the five disease-specific networks. The core network was constructed by the links that appear in four or more of the networks. We mapped the nodes in each network to the family level and then calculated the pan and core of the links. Links with the same source nodes and target nodes and the interaction type (positive or negative) were merged. Links with the same source and target nodes but different interaction types were removed.

To quantify the characteristics of networks, we calculated node attributes such as degree, importance centrality, betweenness centrality, and network attributes such as average path length, modularity index, and scale-free index. The degree, importance centrality, and scale-free index were calculated with an in-house R script with the following method.

We denote d(s_j) as the node degree, defined as the total edges that connect to the node s_j :

We defined the node importance score (importance centrality) to denote the importance of a node by modifying the clustering coefficient (CC). Therefore, the importance score is defined as:

a(s_j ) is the number of edges that exist among the neighbor nodes that directly link to the node s_j , not including s_j itself. b(s_j ) is the maximum number of possible edges among the neighbor nodes of s_j . This definition of CC is as the following:

The reason to modify CC is that the numerator and denominator may be both 0 in φ(s_j ). Moreover, we want CC to be directly proportional to the importance of a node. That is, the more neighbor nodes of s_j and fewer edges among these neighbor nodes, which means the more important of node s_j , the larger the c(s_j ).

Betweenness: the node betweenness is defined by the number of shortest paths going through the node. It is calculated by the betweenness function of the “igraph” package in R (Csardi and Nepusz, 2006).

Scale-free index: if a network has the property where the low-degree nodes are in the majority and hub nodes of high degree are in the minority, the network is defined to be scale-free. The distribution of node degree in a scale-free network follows a power-law distribution, that is,

Thus, the scale-free index is defined as the fitness (r-square) of the linear regression model between log(P(d)) and log(d).

Modularity index: the modularity is calculated by cluster_walktrap and modularity function of the “igraph” package in R (Csardi and Nepusz, 2006), which utilize the random walk algorithm to cluster nodes in the networks.

Average path length: the average path length of the network is calculated by the mean_distance function of the “igraph” package in R.

We first annotated all genes of the microbiome with the KEGG database and the ALCME database separately (Leplae et al., 2010; Wang et al., 2015; Kanehisa et al., 2016). Since the KEGG functional annotation rate for viral genes was low (overall annotation rate: 46.8%), we used annotation of ALCME in later analyses (overall annotation rate: 80.1%). However, the majority of families in ACLAME were without GO or MeGO annotations, and thus we manually annotated them with the protein names that appeared most times within that family. After aggregating the families with the same annotation terms, a total of 2,162 function items were obtained. Thus, we further grouped these items into 16 categories manually by reference to COG (Tatusov et al., 2000). The details of the classification of the 16 categories are described in Supplementary Methods. We further built a database for the detected viral genes with manually classified functional annotation named VirGenFunD. Each sequence was labeled with annotations of taxonomy, KEGG, ACLAME, and VirGenFunD category and classified into five classes.

The enrichment analysis was based on Fisher’s exact test. Here, we use all genes of viruses that take part in a network as the background. Since genes were both annotated with ACLAME family information and taxonomic information, the function of the VirGenFunD category can be retrieved through taxonomic annotation of genes. We first counted the number of each VirGenFunD category in the background list of viruses and then tested each VirGenFunD category for their enrichment in a subset of the background list by counting the number of a VirGenFunD category in that subset and performing Fisher’s exact test with a contingency table (take Category01 as an example):

Numbers that follow the ± sign in the manuscript indicate standard deviations. The Mann–Whitney U test was conducted with the wilcox.test function in R (two-sided), with a significance level set as p ≤ 0.05 for viruses and bacteria. The Kolmogorov–Smirnov test was used to test the normality of the distribution for continuous variables. The t-test was conducted with the t.test function in R (two-sided) for normally distributed continuous variables, with the significance level set as p ≤ 0.05 for metabolites. Fisher’s exact test was conducted with the fisher.test function in R. The networks constructed by SparCC were filtered with r ≥ R, and FDR-adjusted p ≤ 0.05, where r is the SparCC correlation coefficient; R was set as 0.6 in the IBS dataset and 0.4 in the other four datasets. The cutoff is different for the IBS dataset because the sample size is smaller than the other datasets. The network constructed is sensitive to sample size, and to make the network constructed more reliable and get the network of comparable size from different datasets, we set a higher cutoff for the IBS dataset. Statistics for disease-specific network selection were described above.

We used the permutation test to check whether negative or positive correlations were enriched in the relationships between viruses and bacterial phyla. For a given network, we first calculated the number of all bacteria–virus links (n1) and the number of negative bacteria–virus links (n2). For a given bacterial phylum (for example, Firmicutes), we calculated the number of negative Firmicutes–virus links (n3). We sampled n2 links from all the bacteria–virus links (n1) and calculated the number of Firmicutes–virus links (n4) 100,000 times. p-value was then calculated as the upper or lower quantile of n3 in the bootstrapped set of n4. The significance level was set as FDR p ≤ 0.05.

To test the significance of network structures (modularity, scale-free index, and average path length), we built networks of the null models through randomization of preserved number of nodes and edges. The p-values were calculated based on the 1,000 randomized null networks.

The major codes of our analyses are available at https://github.com/lkyvirrrr2001/viruses_analyses/.

The alterations in gut bacteria between IBS and healthy controls have been reported in our previous work (Liu et al., 2016; Wang et al., 2019; Xu et al., 2020). Herein, we further explore the roles of the viral elements using these sequenced metagenomic samples, including 22 IBS patients (diarrhea-predominant) and 15 healthy controls, and 18 paired host metabolomic samples. The complete analysis pipeline is shown in Figure 1A and described in Materials and Methods. Our results demonstrated that the overall ratio of viral gene abundance to all the genes was 9.6% (±3.6%), with no significant difference between cases and controls. A total of 291 viral genera and 50 viral families were detected. Among the annotated viral families, Siphoviridae, Myoviridae, and Podoviridae, all belonging to the Caudovirales order, were the most abundant phages (Figure 1B), which are consistent with the typical composition of the phageome of human adults (Minot et al., 2011; Manrique et al., 2016; Ma et al., 2018). For eukaryotic viruses, Megavirales, Pithoviridae, Baculoviridae, Nudiviridae, Nimaviridae, Circoviridae, Retroviridae, Togaviridae, etc., were detected but were in the minor part. These viruses were also common in the human gut (Table S1) (Mukhopadhya et al., 2019).

The comparison between cases and controls revealed 22 and 6 viral genera depleted and enriched respectively in the IBS group (Mann–Whitney U test, p ≤ 0.05), detected by MetaComp (Zhai et al., 2017). In the previous study, the gut bacteria and host serum metabolites were found to have shifted between IBS and healthy controls (Xu et al., 2020). To explore whether these different viruses were related to the shifts of bacteria and metabolome between IBS and healthy controls, we calculated the Spearman correlations between different viruses and different bacteria (Mann–Whitney U test, p ≤ 0.05) as well as metabolites (t-test, p ≤ 0.05), respectively (Figures 2A, B). IBS-depleted viruses positively correlated with most of the IBS-depleted bacteria, and IBS-enriched viruses positively correlated with IBS-enriched bacteria (Figure 2A), showing a covarying relationship between shifted viruses and bacteria. Among the shifted bacteria, three lactic acid bacteria, Lactobacillus, Lactococcus, and Enterococcus, which are probiotics and beneficial to human health (Hatti-Kaul et al., 2018), were depleted in the IBS group (Mann–Whitney U test, p ≤ 0.05, Table S2). Alongside this shift, C5virus (phage that infects Lactobacillus), Sk1virus (phage that infects Lactococcus), and Phifelvirus (phage that infects Enterococcus) were also depleted in the IBS group (Mann–Whitney U test, p ≤ 0.05) and showed a significant positive correlation with their host bacteria (Spearman correlation, FDR p < 0.01), indicating a lysogenic relationship between the phages and their bacteria hosts.

The host non-targeted metabolomics profiling further supported the relationship between viruses and the disease. A total of 77 negative and 59 positive metabolic ions were detected as significantly different between cases and controls (t-test, p ≤ 0.05). These metabolic ions mostly came from lipids, amino acids, dipeptides, organic acids, disaccharides, benzenoids, etc. (Table S3), which can significantly separate samples of cases and controls by a partial least square discriminant analysis (PLS‐DA) model (Figure S1). Consistency between viruses and metabolome was observed herein, for IBS-depleted viruses positively correlated with IBS-depleted metabolic ions and negatively correlated with IBS-enriched metabolic ions (Figure S2). We also calculated the correlations between viral genes (annotated by ACLAME database) and the different metabolites, and some viral functions significantly correlated with the shifted metabolites (Figure 2B). Some of the correlated viral gene functions deal with viral propagation, such as replication, transcription, and transposition of related functions and structural proteins; others deal with enzymes that are important in biological functions for viruses or their host, such as virally encoded metal-dependent hydrolase, which catalyzes the hydrolysis of a wide range of biologically important substrates including carbohydrates, peptides, and nucleotides. Furthermore, the differential metabolites were involved in 47 metabolic pathways (Table S3), in which 25 overlapped with pathways found in the metagenomic data, and the differential viruses had a significant contribution to these pathways (11 pathways out of 25, Fisher’s exact test, p = 1.67 × 10^-2). These results indicated that the changes of viruses accompanied changes in bacteria and metabolites of the IBS group, and the genes encoded by the viruses may affect the metabolic pathways related to human health.

To get a better understanding of the gut virome in health and diseases, we collected metagenomic datasets from four other diseases, T2D, CD, CRC, and LC, which are typical metabolic or dysfunctional bowel diseases that are marked with dysbiosis in gut microbiota (Qin et al., 2012; Qin et al., 2014; Lewis et al., 2015; Yu et al., 2017). To explore the roles of viruses in a systematic view of diseases, we also collected published metagenomic datasets related to the above four diseases. Meta-data including age and sex of the involved individuals are listed in Table S4. All the raw data from the 438 gut metagenomic samples (including 37 samples in the above IBS cohorts) made up the total 5T size of fastq files which then followed a uniform analysis pipeline (Figure 1A). In the whole gene repertoire of the gut microbiome of all samples, viruses contributed 9.4% (±5.6%) of gene abundance to the total annotated genes (Figure 3A), and no significant difference in the viral ratio was found between each group of cases and controls. The compositions of viruses in the four datasets resembled that in the IBS dataset (Figure 3B), and identifiable phages occupy 68.5%–73.5% of all the viral abundance in different groups (Figure S3). In a systematic perspective of five diseases, there were 106 viral genera (Table S5) and 20 viral families (Table S6) that showed different abundances between cases and controls (Mann–Whitney U test, p ≤ 0.05). However, we found no shared different genus or family among the five datasets.

To explore the possible positive or negative roles of viruses in the human gut with analysis of the disease-related microbiota samples herein, we thus investigated how phages and eukaryotic viruses respectively interacted with bacteria in disease-specific networks. We first constructed the co-abundance network of the microbes at the genus level within each case and control group (see bacterial and viral genus abundance matrixes in Tables S2 and S5). Then, we selected links that specifically appeared in the case group as the disease-specific network (see Materials and Methods). We also constructed healthy-specific networks relative to each case group to make comparisons. Metrics to quantify the properties of the networks were calculated, including mean degree, edge number, importance centrality, betweenness, modularity, scale-free index, and average path length.

All the networks showed frequent correlations within and between bacteria, phages, and eukaryotic viruses, and some of the viral nodes were of high degree and high importance centrality (Figure 4A, Figures S4– S7, and Table S7). The family-level pan (pooled) network and core (shared) network showed that phages, notably the three most abundant Caudovirales phages (Siphoviridae, Myoviridae, and Podoviridae families), took important positions in the network since they had both high degrees and high importance centralities (Figure 4B and Figures S8– S10). The construction of the family-level pan (pooled) network and core (shared) network is described in Materials and Methods.

In all healthy-specific interaction networks, the mean degrees of viruses were comparable to those of bacteria. The mean degrees of eukaryotic viruses were even higher than those of bacteria in the healthy control group of the T2D dataset (Table S8, Mann–Whitney U test, W = 8079.5, p = 0.01), while in T2D, CD, and LC, the mean degrees of phages or eukaryotic viruses were smaller than those of bacteria which might indicate the decreased number of relationships of viruses with other microbes in these disease-specific networks (Table S8).

We then focused on the relationships between viruses and bacteria in disease-specific interaction networks to find how viruses impacted different bacteria phyla. Most of the relationships linked to viruses were from bacteria, indicating that viruses and bacteria had close and intricate interactions (Figure 4A and Table S9). Specifically, viruses, including both phages and eukaryotic viruses, had rarely negative relationships within themselves, compared to a high negative correlation ratio between viruses and bacteria (Figure 4C and Table S9, proportion test, FDR p < 0.05). This was also the case when we combined all the groups into a pan-group (Figure S11). The ratios of negative correlations in phage–phage correlation, phage–eukaryotic virus correlation, and eukaryotic virus–eukaryotic virus correlation (14.0%, 5.5%, and 3.1%, respectively) are significantly lower than the ratios of negative correlations in bacteria–phage correlation, bacteria–eukaryotic virus correlation, and bacteria–bacteria correlation (37.1%, 42.4%, and 25.9%, respectively) (proportion test, FDR p < 0.05). This result suggested that viruses might have restraints on some bacteria and rarely conflicted within themselves. Firmicutes, Bacteroidetes, Proteobacteria, and Actinobacteria were four major bacterial phyla that interacted with viruses, while Firmicutes and Bacteroidetes had the highest interaction ratio and thus might be the most impacted by viruses (Figure 4D).

Moreover, viruses showed different preferences in positive or negative correlations with these four bacteria. Viruses tended to have more positive correlations with Firmicutes and Actinobacteria and more negative correlations with Proteobacteria and Bacteroidetes (permutation test, Figure 4E), which are common dysbiosis-related and inflammation-mediating bacteria, respectively (Shin et al., 2015; Wassenaar and Zimmermann, 2018). This phenomenon was even more significant in healthy-specific networks (Figure S12), indicating a possible positive role of viruses to inhibit disease-mediating bacteria in both health and diseases.

To further characterize the roles of viruses in the interaction networks of the disease-related microbiome, we calculated network structural features including modularity and scale-freeness of networks, to explore how viruses impacted the network structure. Modularity is a measurement of the property that a network can be divided into individual communities in which members are densely interconnected and sparsely connected outside. The modularity of biological networks makes the subgroups function semi-autonomously (https://psychology.wikia.org/wiki/Modularity_(biology)). In all healthy-specific and disease-specific networks, the modularity indices were significantly higher than those of random network null models with the same number of nodes and edges of each network (1,000 times of randomization, p < 0.001) (Figure 5A and Figures S13– S16). These differences showed that relative independent communities existed among microbes, in which viruses interspersed among bacteria rather than grouped into individual modules. More interestingly, we explored the viral effects on the structures of networks by removing the virus nodes in disease-specific networks and observed that the modularity of the bacteria network increased in all groups (Figure 5B and Table S9), meaning that viruses connect between bacteria modules and contribute to shaping the network structure of the gut microbiota.

The scale-free property of networks, in which a few nodes possess a large number of relationships, and most nodes possess a small number of relationships, is typical of biological systems that are robust to random disruptions (Barabási and Albert, 1999). Herein we calculated scale-free indices, defined as the r-square of the node degree fitting a log-transformed power-law distribution, of the five disease-specific interaction networks (Table S9). All the networks showed the characteristics that low-degree nodes had higher frequencies than high-degree nodes (Figure 5C, 1,000 times of randomization, p < 0.001), and short paths existed among different taxa (Table S9), compared to the average path length of six in random networks (1,000 times of randomization, p < 0.001) (Watts and Strogatz, 1998; Faust et al., 2012). This result indicated that the disease-specific networks had the property of scale-freeness. When comparing the scale-free indices of the total networks with the ones in which viruses were not considered (non-viral), we found that the scale-free indices of all the networks decreased (Table S9, 1,000 times of randomization, p < 0.001). Although the decrease in scale-free indices was mainly attributed to the removal of substantial amounts of low-degree viral nodes, this result indicated that viruses improved the scale-freeness of the networks and made the microbial community of gut microbiota more robust.

The network structure analysis illustrated that viruses could bridge between bacteria, but not all viruses were equally important. To determine which viruses were more important in communicating between bacteria, we defined the ones that could shorten all the bacteria’s average path lengths as key viruses. By calculating the betweenness centrality of each viral node and adding them to the bacterial community in decreasing order, we found lists of such key viruses that minimize the average path length in the pan network as well as five disease-specific networks (Figure 5D and Table S10). These viruses took significant parts in interacting with short-chain fatty acid-producing bacteria Firmicutes (121 out of 125 interacting bacteria, Fisher’s exact test, p = 1.58 × 10^-2) and Bacteroidetes (79 out of 83 interacting bacteria, Fisher’s exact test, p = 6.23 × 10^-2), which implied that in the disease-specific networks, viruses might take an important part in the gut metabolism, by the coaction with these bacteria.

To further explore how viruses interact with bacteria, we first investigated the viral gene functions by annotating them with viral protein families in the ACLAME database. Since only a small part of protein families are annotated with GO or MeGo, we reannotated them and categorized them into 16 categories (Figure 6 and Table S11). Based on the 16 categories, we performed enrichment tests on gene functions of viruses that interact with bacteria (see Materials and Methods). The gene set of key viruses selected in each network described above was generally enriched in functions of “Transporter activity”, which mainly includes ATP-binding cassette (ABC) transporter- and transmembrane transporter-related proteins, “Chaperons and secretion system”, which mainly includes proteins involved in type III and type IV secretion system that transport molecules from bacterial cells to other cells, “Metabolic enzymes”, which include general enzymes such as oxidoreductase, hydrolase, and modification-related enzymes, and “Signal transduction”, which mainly includes a two-component signal transduction system and response to stress (Table S12, Fisher’s exact test, FDR p < 0.01). These functions are essential biological functions for both viruses and hosts. Viruses in these ways might affect the phenotype of bacteria under the condition of diseases.

We further tested the significance of functional enrichment of phages that negatively and positively interacted with bacteria (phages that had negative link ratios of 1 and 0 were selected as representatives, respectively). In comparison with all the annotated phage genes, functions that were typical of temperate phages such as “Integration and recombination” and “Metabolic enzymes” were enriched in both positive and negative interacting phages with bacteria (Tables S13 and S14, Fisher’s exact test, FDR p < 0.01), supporting that most phages in the gut were temperate phages (Toussaint and Rice, 2017). Function “Transcription and regulation of gene expression” and “Nucleotide transport and metabolism” were only enriched in phages that positively interact with bacteria. In contrast, function “Toxins and detoxification” (this category includes toxins, post-segregating killing process, and detoxification of mercury ion) was only enriched in phages that negatively interact with bacteria. These results indicated the possible mechanisms underlying the correlations between phages and bacteria, either through direct predation or lysogeny, or through indirect interaction between phages and non-host bacteria by affecting the host gene expression, metabolism, and virulence. We performed the same analysis in the healthy-specific networks but found no positive-interaction-specific or negative-interaction-specific viral function, indicating a homogeneous distribution of gene functions in positive and negative interacting phages with bacteria (Tables S15 and S16). As for eukaryotic viruses, which were generally ignored in the gut virome, they interacted with bacteria with the enriched functions of “Signal transduction”, “Metabolic enzymes”, and “Transporter activity” (Tables S17 and S18, Fisher’s exact test, FDR p < 0.01) in both disease- and healthy-specific networks, giving hints to the potential roles of eukaryotic viruses in the gut.

To provide convenience for further analyses and more insights into the role of viruses in the human gut, we summarized our detected viral sequences and annotations from different sources into a database named VirGenFunD. Although there are some excellent viral databases like GVD (Gregory et al., 2020), ACLAME (Leplae et al., 2004), and pVOGs (Grazziotin et al., 2016), which include large amounts of viral sequences, they either lack functional annotations or have sporadic GO annotations. The VirGenFunD database contained 3,351,765 viral gene sequences from five disease-related datasets and the corresponding healthy controls. The VirGenFunD annotation thus nearly doubled the number of the known function categories (Table 1) and provided vital clues for understanding the functions of viruses in the human gut. The VirGenFunD database is available at http://cqb.pku.edu.cn/ZhuLab/VirGenFunD/, https://yjiang724.github.io/VirGenFunD/.