Published: February 12, 2020

Introduction {#sec1}
============

Malnutrition is a major health concern in low- and middle-income countries and the leading cause of death in children younger than 5 years ([@bib46]). Child malnutrition falls under three main subtypes: underweight, stunting, and wasting. Stunting is defined by the World Health Organization (WHO) as a height-for-age *Z* score (HAZ) that is two standard deviations below the WHO Child Growth Standards mean or more ([@bib30]). The onset of stunting occurs in children between the ages of 6 and 23 months old, when linear growth is highly susceptible to nutritional deficiency and environmental stress. This is also the transition period from exclusive milk feed to solid food. Most children from countries with high stunting rates do not have regular access to nutrient-rich or nutrient-balanced foods, which emphasizes the role of diet in child stunting ([@bib37]). In addition to nutritional deficiencies, stunting is linked to repeated diarrheal infections and poor sanitation ([@bib33]). This multifactorial condition has multigenerational consequences and is associated with increased susceptibility to chronic diseases, poor education performance, and poor socioeconomic conditions ([@bib4], [@bib68]). Although the health risks and environmental drivers of child stunting are well documented, the underlying pathophysiological mechanisms for this disease remain largely unknown, and recent studies have unveiled a correlation between an altered gut microbiota and stunting ([@bib39], [@bib68]).

The gut microbiota refers to the community of microorganisms living in the gastrointestinal tract, including bacteria and bacteriophages (or phages). The gut microbiota plays essential roles in host metabolism, immune modulation, and colonization resistance to pathogens. In healthy individuals, approximately 90% of the gut bacterial community is dominated by members of the Firmicutes and Bacteroides phyla, whereas members of the Proteobacteria, containing many pathogens, reside at low levels ([@bib4], [@bib68]). Shifts in gut bacterial composition have been associated with increasing numbers of diseases including inflammatory bowel diseases (IBD), allergies, diabetes, and obesity ([@bib4], [@bib42]). Generally, an increase in members of the Proteobacteria and lower overall community diversity are observed, in addition to disease-specific shifts of the gut microbiota. Stunted children also harbor an altered gut microbiota, with increased prevalence of *Enterobacteriaceae*, which may further contribute to their impaired growth and nutritional deficiencies. Indeed, the local inflammation caused by overgrowth of *Enterobacteriaceae* has been shown to lead to impaired digestive and absorptive functions of the gut, all linked to stunting ([@bib26], [@bib68]). In addition, recent studies further report that stunted children have an immature gut microbiota relative to their non-stunted age-matched counterparts. This immaturity was defined as a lower microbiota-for-age *Z* score, characterized by a lower α-diversity of the gut microbiota and disproportionately higher levels of Proteobacteria ([@bib62]).

In contrast to the work done to date on gut bacterial communities, the role of gut phage communities in child stunting remains largely unexplored. Phages are bacterial viruses that are key to the maintenance and function of many ecosystems by supplying bacteria with genes involved in host adaptation, toxin production, and metabolism ([@bib59], [@bib65]). Phages control bacterial diversity and abundance and can modify the O-antigen component of lipopolysaccharide (LPS) in gram-negative bacteria, which is particularly relevant to intestinal inflammation ([@bib67]). In the gut, phages are abundant, with a phage-to-bacteria ratio close to 1:1 based on sequencing data ([@bib60]). Phage communities are distinct between individuals and stable over time compared with gut bacterial ones, despite a few examples of common and widespread phages, such as crAssphage, and rapidly evolving phages ([@bib40], [@bib69]). Recently, changes in the diversity and abundance of phages have been associated with many diseases such as IBD, diabetes, malnutrition, AIDS, and Parkinson's disease, highlighting their potential role in human health ([@bib16], [@bib31], [@bib41], [@bib48], [@bib64]).

Despite the known effects of phages on bacterial communities, bacteria-phage dynamics in the gut remains poorly described ([@bib18], [@bib35], [@bib60]). This is an important limitation to our functional understanding of how the gut microbiota alters human health. Relevant examples of the narrow host range of phages and their effective therapeutic use provide compelling arguments that phages could be used to manipulate gut bacterial communities to promote health ([@bib24], [@bib59]). Here, we explored whole community bacteria-phage population dynamics from stunted and non-stunted, otherwise clinically healthy, children from Dhaka, Bangladesh. Despite the significant progress in tackling malnutrition, Bangladesh still has one of the highest rates of child stunting in the world, affecting an average of 36.1% of children younger than 5 years ([@bib47]). We obtained fecal samples from 30 non-stunted and 30 stunted children from the Mirpur slum area, aged between 14 and 38 months, and separated the free phage and bacterial communities. Because diet is a major regulator of gut bacterial diversity, we split the children into two age groups based on their age-specific diet, with each group consisting of 15 non-stunted and 15 stunted children. Children in the younger group, aged between 14 and 23 months, had a mixed diet of milk and solid food, whereas the older children, aged between 28 and 38 months, had fully transitioned to solid food for more than 2 months. After characterizing the abundance and diversity of the phage and bacterial communities in the original fecal samples, we proceeded with *in vitro* 2 × 2 factorial cross-infections of both communities, where bacteria from non-stunted children were exposed to phages from stunted children and vice versa. This allowed us to (1) identify the differences in bacteria and phage communities between non-stunted and stunted children; (2) determine how phages modify gut bacterial communities of non-stunted and stunted children; and (3) investigate the role of the children's age, diet, and health conditions on bacteria-phage interactions.

Results {#sec2}
=======

Children Metadata {#sec2.1}
-----------------

The younger age group of children comprised an equal number of boys and girls with a mean age of 19.5 ± 2.5 months, mean weight of 9.6 ± 1.3 kg, mean height of 77.7 ± 3.4 cm, and a mean HAZ score of −1.8 ± 0.9. Children in the older age group comprised 18 girls (60%) and 12 boys (40%), with a mean age of 33.2 ± 2.9 months, mean weight of 11.5 ± 0.9 kg, mean height of 87.4 ± 3.1 cm, and a mean HAZ of −1.7 ± 0.92 ([Table S1](#mmc2){ref-type="supplementary-material"}). In younger children, 90% of them were fed breast milk; the rest received a mixture of breast milk and formula. Median HAZ scores of non-stunted children in the younger group (−0.97; IQR: −1.44 to −0.73) and the older age group (−1.27; IQR: −1.47 to −0.55) were significantly higher (Wilcoxon, matched-pairs signed rank test, p \< 0.0001) than the median HAZ scores of the children with stunting in the younger age group (−2.53; IQR: −2.97 to −2.27) and the older stunted children (−2.32; IQR: −2.79 to −2.1), respectively.

Bacterial and Phage Abundances in Non-stunted and Stunted Children {#sec2.2}
------------------------------------------------------------------

To examine differences in bacterial and phage absolute abundances among all health status and age groups, we counted stained bacteria and viral-like particles (VLPs) by microscopy (see [STAR Methods](#sec4){ref-type="sec"}) ([Figure S1](#mmc1){ref-type="supplementary-material"}A). There were no significant differences in bacterial or VLP abundances, irrespective of health status or age group ([Figures S2](#mmc1){ref-type="supplementary-material"}A and S2B). Viral abundances in young children were more variable compared with the older ones, regardless of their health status ([Figure S2](#mmc1){ref-type="supplementary-material"}A). This variability in VLP abundance is specific to young children and was not observed with the bacterial abundance data. We then determined the resulting virus-to-bacteria ratio (VBR) in all children. We report no statistical difference between the VBRs of the different groups, most probably because of the high variability in VLP abundance in the younger children ([Figure S2](#mmc1){ref-type="supplementary-material"}C).

Bacterial Diversity in Non-stunted and Stunted Children {#sec2.3}
-------------------------------------------------------

We used 16S rRNA gene and shotgun sequencing ([Figure S1](#mmc1){ref-type="supplementary-material"}B; [Table S3](#mmc4){ref-type="supplementary-material"}) to examine how bacterial community composition may differ according to health status and age group. We further examined how these changes in phylum or species composition affected the overall diversity by determining the Shannon Evenness Index, as well as bacterial richness of each age group. Actinobacteria were the most abundant in younger children regardless of health status, though their relative abundance was slightly higher in non-stunted children ([Figure 1](#fig1){ref-type="fig"}A). The second most common phylum in younger children was the Firmicutes. In contrast, in older children, the Firmicutes were the most abundant, followed by the Actinobacteria ([Figure 1](#fig1){ref-type="fig"}A). Stunted children carried more Proteobacteria compared with their non-stunted counterparts, with younger stunted children showing the highest proportion of Proteobacteria (9%; [Figure 1](#fig1){ref-type="fig"}A). At the species level, non-stunted younger children had more *Bifidobacterium adolescentis* (6%), whereas stunted children had a slightly greater proportion of *Bifidobacterium longum* (25%) and notably greater proportions of *Escherichia coli* (19%) ([Figure 1](#fig1){ref-type="fig"}B). For older non-stunted children, we observed a greater proportion of *B. adolescentis* (increasing from 6% to 13%), as well as an increase in *Lactobacillus salivarius* (6%); stunted older children carried more *Klebsiella pneumoniae* (19%) and *E. coli* (16%) compared with their non-stunted counterparts. In general, the number of reads assigned to *B. longum* was greater in younger children than older ones, whereas *Bifidobacterium kashiwanohense* was the dominant species in all age and health groups ([Figure 1](#fig1){ref-type="fig"}B).Figure 1Bacterial Communities Are Less Diversified in Younger, Stunted ChildrenRelative abundance of bacterial phyla based on 16S analysis (A) and species based on shotgun metagenomics (B) in children's fecal samples (n = 60). Samples were pooled for sequencing (see [STAR Methods](#sec4){ref-type="sec"}). Children in the younger age group were aged between 14 and 23 months; children in the older age group were aged between 23 and 38 months. Shannon Diversity Index was compared among (C) age groups according to health status, and (D) HAZ scores according to age groups.See also [Figures S2](#mmc1){ref-type="supplementary-material"} and [S3](#mmc1){ref-type="supplementary-material"}.

Though our sample size was too small to test for significance, we noted that the difference in bacterial species diversity and richness between non-stunted and stunted children tended to be greater for the older age group than the younger one ([Figures 1](#fig1){ref-type="fig"}C and 1D). This was due to a greater increase in diversity in non-stunted children with age. At the phyla level, diversity and richness only increased in the non-stunted group ([Figure S3](#mmc1){ref-type="supplementary-material"}). Overall, this suggests that the effects of stunting became more apparent with age ([@bib4], [@bib15], [@bib68]).

Abundance of Proteobacteria Pathogenic Species {#sec2.4}
----------------------------------------------

Given the suggested link between pathogenic bacteria in the gut and child stunting ([@bib20], [@bib68]), we analyzed the relative abundance of known human pathogenic species, focusing on Proteobacteria species (see [STAR Methods](#sec4){ref-type="sec"}). Within the Proteobacteria, *E. coli* was the most abundant potential pathogenic bacterial species in all age and health groups, except in older stunted children, where *Klebsiella pneumonia* dominated at 51% of the pathogenic community, followed by *E. coli* at 44% ([Figure 2](#fig2){ref-type="fig"}A). Young stunted children had the highest occurrence of *E. coli* (93% of reads within Proteobacteria) overall. We further analyzed the presence of pathogenic and non-pathogenic *E. coli* strains based on their nucleotide identity to known *E. coli* pathogenic strains ([Table S2](#mmc3){ref-type="supplementary-material"}). The ratio of pathogenic to non-pathogenic strains of *E. coli* was higher for younger children ([Figure 2](#fig2){ref-type="fig"}B). Older stunted children had a relatively lower abundance of *E. coli*, likely due to competition with *K. pneumonia*, which dominated this age group ([@bib22]). *E. coli* BW25113 and *E. coli* ER2796 were the dominant strains in older stunted children (Kruskal-Wallis, Dunn's post hoc test, p \< 0.05), whereas *E. coli* LY180 and *E. coli* O104:H21 strain CFSAN002236 were the major *E. coli* strains in their non-stunted counterparts (one-way analysis of variance \[ANOVA\], Holm-Sidak's post hoc test, p \< 0.05) ([Figures 2](#fig2){ref-type="fig"}C and 2D).Figure 2Non-stunted and Stunted Children Harbor Many Pathogenic *E. coli* Strains(A) Relative abundance of known pathogenic species of Proteobacteria.(B) Relative abundance of pathogenic and non-pathogenic *E. coli* strains.(C and D) *Z* scores of the abundance of *E. coli* strains in children from the younger (C) and older (D) age groups. When present, error bars correspond to the standard deviation of bacterial abundances.

Phage Diversity in Non-stunted and Stunted Children {#sec2.5}
---------------------------------------------------

The majority of reads we identified as viruses were phages (94% on average); reads identified as nonbacterial viruses varied from 0.5% for younger stunted children to 14% for older non-stunted children. Nonbacterial viruses were archaeal, human, and environmental viruses including amoeba and algae viruses. The proportion of nonbacterial viruses was highest for non-stunted older children, with most reads belonging to Circoviruses, a DNA virus commonly found in humans ([@bib29]). On average, 36% of the reads could not be assigned to any reference phage genomes.

Members of the Siphoviridae family dominated phage reads in all children ([Figure 3](#fig3){ref-type="fig"}A). Young, non-stunted children had high proportions of Podoviridae (28%), Myoviridae (14%), and Microviridae (6%); whereas older, non-stunted children had a greater proportion of reads assigned to Myoviridae (21%) and Microviridae (18%). In the stunted group, young children had less Podoviridae (12%) than older ones (23%). At the species level, *Bacillus* phages were the most abundant phages in young, non-stunted children, totaling 14% of the reads ([Figure 3](#fig3){ref-type="fig"}B), followed by crAssphage at 7%. Young, stunted children carried more *Lactococcus* phages (56%), followed by *E. coli* phages (6%). In contrast, Microviridae and *Bacteroides* phages were the most abundant phages in older, non-stunted children (both at 13%), closely followed by *Bacillus* phages (12%), whereas older, stunted children had more *Salmonella* (21%) and *Klebsiella* (13%) phages. In terms of diversity and richness, phage species followed a pattern similar to that of bacteria, with a tendency for higher diversity in non-stunted children compared with stunted ones ([Figures 3](#fig3){ref-type="fig"}D and 3E). Contrasting with the bacterial data, the difference in phage diversity between non-stunted and stunted children tended to decrease with age ([Figure 3](#fig3){ref-type="fig"}D).Figure 3Phage Community Composition and Diversity Indices Change with Age and HealthRelative abundance of phage species in all fecal samples. Samples were pooled for sequencing (see [STAR Methods](#sec4){ref-type="sec"}). Relative abundance of (A) phage families, (B) phage species, and (C) abundance of temperate phages presented as *Z* score. Shannon diversity in fecal samples of children was compared among (D) age groups according to health status, and (E) HAZ scores according to age groups.See also [Figure S2](#mmc1){ref-type="supplementary-material"}.

We next aimed to determine the dominant phage replication cycle, as the effects of phages can vary based on their replication cycle (lytic versus lysogenic cycle) ([@bib42]). Most temperate phages integrate their genome into their bacterial host, thus protecting their host against further infection by closely related phages, and potentially transferring genes providing a competitive advantage ([@bib60]). We identified temperate phages based on the presence of the *integrase* or *integrase*-like genes, which are used by most known temperate phages (see [STAR Methods](#sec4){ref-type="sec"}). On average, young children had more temperate phages than older children, irrespective of health status; however, within each age group, non-stunted children had more temperate phages than their stunted counterparts (young: non-stunted *Z* score = 1.2 and stunted *Z* = −1; old: non-stunted *Z* = 0.4 and stunted *Z* = −0.5) ([Figure 3](#fig3){ref-type="fig"}C).

Bacterial and VLP Abundances and Diversity after Cross-Infections {#sec2.6}
-----------------------------------------------------------------

To gain a whole community perspective on bacteria-phage interactions, and to test if gut phages isolated from non-stunted and stunted children alter gut bacterial communities from children with different health statuses, we proceeded with fully factorial cross-infection experiments in each age group ([Figure S1](#mmc1){ref-type="supplementary-material"}). We then determined bacterial and viral absolute abundances ([Figure 4](#fig4){ref-type="fig"}), as well as bacterial and phage community composition. While the VBRs in younger children remained similar between the start and the end of the cross-infections ([Figure 4](#fig4){ref-type="fig"}A), we noted a significant decrease in VBRs in older children between the original fecal samples and the cross-infected cohorts ([Figure 4](#fig4){ref-type="fig"}D; p \< 0.05, one-way ANOVA, Holm-Sidak's multiple comparisons test). This change was mostly driven by an increase in bacterial abundances during incubation ([Figures 4](#fig4){ref-type="fig"}E and 4F).Figure 4Virus-to-Bacteria Ratios Change after Cross-Infections in an Age-Specific MannerViolin plots of whole community bacterial and VLP abundances in children from the younger (A--C) and older (D--F) age groups prior and post cross-infections. (A and D) Virus-to-bacteria ratios; (B and E) bacterial abundances by epifluorescence microscopy; (C and F) abundance of VLPs by epifluorescence microscopy. Bars connected by the same letter are not significantly different (p \< 0.05, one-way ANOVA, Holm-Sidak's multiple comparisons test). ND, not detected; N, non-stunted; S, stunted; b, bacteria; p, phage; HK, heat killed. The width of the plots reflects the data frequency distribution; dotted lines show the median and lower and upper quartiles.

The larger sample size for the cross-infections allowed us to test for significant differences in diversity and richness among cohorts. We found that at the phyla level, diversity was significantly greater in cross-infected cohorts of younger children relative to older children ([Figure S4](#mmc1){ref-type="supplementary-material"}A; t = 7.4, p \< 0.00001). Bacterial phyla richness was greater in older stunted children than older non-stunted children (t = −7.78, p = 0.001) after cross-infections ([Figure S4](#mmc1){ref-type="supplementary-material"}B), whereas younger stunted and non-stunted children had similar richness (t = 0, p = 1). When bacteria from non-stunted children (Nb) were exposed to phages from stunted children (Sp), Proteobacteria (mainly *Enterobacteriaceae*) increased by 6% relative to the heat-killed control (Nb-HKSp) and nearly doubled from the incubation control (NbNp) in younger children (treatment: 6,719.6 ± 8,818.7 reads; heat-killed control: 3,847 ± 3,659.1 reads; incubation control: 3,223.6 ± 2,827.9 reads) ([Figure 5](#fig5){ref-type="fig"}A). In contrast, Proteobacteria decreased when bacteria from stunted children (Sb) were exposed to phages from non-stunted children (Np). This did not seem to be driven by the phages themselves, as the difference between live phages or heat-killed phages was small (treatment: 2,279 ± 1,609.7 reads; heat-killed: 1,621.6 ± 2,315.3 reads; incubation control: 6,223 ± 8,054 reads). A similar trend was observed for Firmicutes (mostly *Veillonellaceae*) when bacteria from stunted children were exposed to phages from non-stunted children: the relative abundance of Firmicutes increased compared with that of the controls (treatment 11,650.6 ± 8,396.1 reads; heat-killed: 5,688 ± 3,680.4 reads; incubation control: 6,546 ± 3,386.9 reads). Such changes in dominant phyla were not observed in older children.Figure 5Phage Communities Alter Bacterial Diversity in an Age-Specific Manner(A) Relative abundance of bacterial phyla after cross-infections in both age groups(B) Bacterial species contribution to β-diversity (SCBD) in younger (left) and older (right) children. Shown are the top 30 species with greatest SCBD values for each age group, where green points correspond to species with SCBD values greater than the mean SCBD values of all species (n = 179 species total).(C) Relationship between local contribution of samples to bacterial β-diversity (LCBD) and diversity indices after cross-infections for the younger (top) and the older age groups (bottom). Graphs 1 and 4 correspond to LCBD versus treatment type, where cross-infection samples with non-stunted bacteria are shown in red, and those with stunted bacteria are in green. Graphs 2 and 5 correspond to LCBD versus species richness. Graphs 3 and 6 correspond to LCBD versus Shannon evenness.See also [Figures S4](#mmc1){ref-type="supplementary-material"} and [S5](#mmc1){ref-type="supplementary-material"}.

Bacteria and Phage Species Turnover after Cross-Infection {#sec2.7}
---------------------------------------------------------

To identify which species varied the most among the cross-infection treatments, we quantified the contributions of each species to β-diversity (or species turnover; calculated here as the species contribution to β-diversity \[SCBD\]) for both the bacteria and phage species. Similarly, to test which treatment had the most unique community composition and greatest turnover relative to other treatments, we calculated the contribution of each treatment to the total β-diversity (or local contribution to β-diversity \[LCBD\]). For the bacterial species of younger children, the only treatment with a significant LCBD was bacteria from non-stunted children (Nb) exposed to phages from stunted children (Sp) relative to their controls. This treatment also had the lowest bacterial Shannon Diversity Index compared with its controls ([Figure 5](#fig5){ref-type="fig"}C). The species that contributed the most to change was *Amycolatopsis lurida* (highest species contribution to β-diversity; [Figure 5](#fig5){ref-type="fig"}B). When bacteria from young stunted children (Sb) were infected with phages from young non-stunted children (Np), we found that the heat-killed controls had a significant LCBD and highest Shannon Diversity, but the lowest richness (thus they were unique because of a loss in species dominance; [Figure 5](#fig5){ref-type="fig"}C). *E. coli* showed the highest SCBD score, followed by *B. longum* ([Figure 5](#fig5){ref-type="fig"}B).

For the phage species of younger children, similar to the bacteria, the treatment with a significant LCBD was the bacterial communities from non-stunted children exposed to phages from stunted children ([Figure S5](#mmc1){ref-type="supplementary-material"}). Contrasting with our bacterial results, this treatment did not have the lowest evenness among the treatment groups. In terms of species contributing the most to this change, we found that the Gokushovirinae, a ssDNA phage associated mostly with Gammaproteobacteria ([@bib57]), had the highest contribution to β-diversity, followed by *Parabacteroides* phages.

The changes in the older children were limited and did not follow those of young children. For the bacterial species, there were no treatments with a significant LCBD ([Figure 5](#fig5){ref-type="fig"}C). However, for the phage communities, both treatments showed a significant change (LCBD) relative to their respective controls ([Figure S5](#mmc1){ref-type="supplementary-material"}). The *Sphingomonas* phages contributed the most to species turnover when bacteria from non-stunted children were infected by phages from stunted children, whereas *Lactococcus* phages had the highest score when bacteria from stunted children were infected with phages from non-stunted children. Furthermore, the pathogenic-nonpathogenic *E. coli* ratio did not change after cross-infections specific to different treatments within age groups (not shown).

Multiple Factor Analysis of the Interactions between Phages, Bacteria, and the Metadata {#sec2.8}
---------------------------------------------------------------------------------------

To synthesize the cross-infection results, we examined the correlations among the dominant bacterial and phage species, and the metadata (diet, age, sex, health status, cross-infection, HAZ, WAZ \[weight-for-age Z\], and WHZ \[weight-for-height Z\] of children) using a multiple factor analysis (MFA) framework ([@bib5]). This multivariate correlation allowed us to assess the contribution of each metadata variable to changes in the bacterial and phage populations. The first dimension of the MFA explained 36.37% of the total variation in the phage and bacterial community and descriptor variables ([Figure 6](#fig6){ref-type="fig"}A), and separated variables according to their correlation with age ([Figure 6](#fig6){ref-type="fig"}B). The second dimension of the MFA explained an additional 21.13% of the total variability in all three data matrices and largely separated variables correlated to non-stunted versus stunted children ([Figure 6](#fig6){ref-type="fig"}B). In light of these two axes (axis 1: age gradient, axis 2: health status gradient), we found that among the dominant bacteria and phage species, certain species tended to be more prominent in older children (i.e., *Salmonella* phage and *Bacteroides dorei*), whereas others were more prominent in the younger children (i.e., Gokushovirinae phage and *E. coli* and *B. longum*; [Figure 6](#fig6){ref-type="fig"}A). The Salmonella phage was also associated with the temperate phage group. Finally, certain phage taxa (*Lactococcus* and *Parabacteroides*) were more common in non-stunted children, whereas Gokushovirinae were more prominent in stunted children after cross-infections.Figure 6Multiple Factor Analysis of the Dominant Bacterial Species, Phage Species, Phage Replication Cycle, and Children Metadata for Both Age GroupsDominant bacterial and phage species shown were selected using a PCA and the contribution circle from the cleanplot.pca() function.(A) MFA ordination of quantitative variables, with the contribution of each qualitative variable to dimensions 1 and 2 of the MFA as insets.(B) MFA ordination of qualitative variables. Quantitative variables include milk feed, age, sex, HAZ, WAZ, and WHZ; qualitative variables include health status, treatment type, and treatment factor. The Hellinger transformation of the bacterial and phage data, and the standardized values for the temperate replication cycle are presented. Quantitative children metadata was scaled.See also [Figures S6](#mmc1){ref-type="supplementary-material"} and [S7](#mmc1){ref-type="supplementary-material"}.

We also ran MFAs to examine changes in phage and bacterial communities within each age group ([Figures S6](#mmc1){ref-type="supplementary-material"} and [S7](#mmc1){ref-type="supplementary-material"}), this time examining all phage and bacterial species. In each case, the first two dimensions of the MFAs explained a large portion of the total variance in both communities and descriptor variables (\>60%). In young children, the first dimension separated the non-stunted children with predominantly *Parabacteroides*, *Clostridium*, *Enterobacteria*, *Burkholderia*, and *Lactococcus* phages, and *E. coli*, from the stunted ones who had greater proportions of Gokushovirinae, *Chlamydia* phages, and numerous *Bifidobacterium* ([Figure S6](#mmc1){ref-type="supplementary-material"}). The second dimension of the MFA separated heat-killed controls from treatments. Heat-killed controls were more strongly associated with temperate phages. In older children ([Figure S7](#mmc1){ref-type="supplementary-material"}), the first dimension once again separated the non-stunted from stunted children, and the second dimension separated heat-killed from treated samples. As with the younger age group, the heat-killed control was associated with temperate phages. We also found that the bacterial and phage species associated either with non-stunted (*Bacteroides dorei*, *K. pneumoniae*, and the *Rhodococcus*, *Bordetella*, *Klebsiella*, and *Salmonella* phages) or stunted children differed from those identified in younger children, as noted in [Figure 6](#fig6){ref-type="fig"}.

MFA of the Interactions between Functional Traits and Metadata {#sec2.9}
--------------------------------------------------------------

Our last set of analyses examined changes in functional traits within original samples and within cross-infected cohorts. In line with previous reports ([@bib41], [@bib44]), most metagenome reads (70% on average) mapped to genes with unknown function. In the original stool samples, we determined which Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway was most abundant in each age and health group for both phages and bacteria. To evaluate this, we conducted an MFA correlating the functional trait matrix (expressed as relative abundances to normalize for the total number of reads per sample) with metadata. We only examined the dominant functional traits and removed the unknown traits to simplify the interpretation of the MFA ([Figure 7](#fig7){ref-type="fig"}). As before, we found that at the onset (original stool samples), non-stunted children had higher HAZ and WAZ (weight-per-age) scores than their stunted counterparts, whereas only young, non-stunted children had higher WHZ (weight-per-height) scores relative to other age and health groups ([Figures 7](#fig7){ref-type="fig"} and [S8](#mmc1){ref-type="supplementary-material"}). The MFA showed that young, stunted children had a greater number of traits related to human disease and metabolism, as well as traits related to cell growth and/or death and genetic replication and/or repair compared with non-stunted, older children ([Figure 7](#fig7){ref-type="fig"}). We note however that other (undetermined) pathways related to genetic replication and/or repair were likewise correlated to older children, which in turn were correlated to more solid diet and females. For the cross-infection samples, we were interested in identifying more specific changes in KEGG pathway abundances across treatments relative to their controls. We thus examined the turnover in pathway composition among cohorts, that is, the contribution of treatments (LCBD) and pathways (SCBD) to β-diversity. For the latter, we arbitrarily present the top 30 traits with highest contribution to β-diversity ([Figure 7](#fig7){ref-type="fig"}). In general, traits were most variable among treatments in younger children ([Figure S8](#mmc1){ref-type="supplementary-material"}). Traits that contributed most to β-diversity in the younger age group were related to human disease, cellular process, and membrane transport, whereas those contributing most to β-diversity in the older age group were cellular process, signal transduction, carbohydrate metabolism, and xenobiotics biodegradation and/or metabolism ([Figure S8](#mmc1){ref-type="supplementary-material"}).Figure 7MFA of the Interactions between Bacterial Functional Traits and Metadata(A) MFA ordination of dominant quantitative variables, with the contribution of each quantitative variable to dimensions 1 and 2 of the MFA as insets. Quantitative variables are as in [Figure 6](#fig6){ref-type="fig"}. Functional traits are expressed as relative abundances to normalize for the total number of reads per sample.(B) Functional trait contribution to β-diversity (SCBD), based on KEGG metabolic pathway abundances, for younger (left) and older (right) children. Shown are the top 30 traits with greatest SCBD values for each age group, where green points correspond to traits with SCBD values that are greater than the mean SCBD values of all traits (n = 8 and 6 for younger and older children, respectively).See also [Figure S8](#mmc1){ref-type="supplementary-material"}.

Discussion {#sec3}
==========

Child stunting affects 22% of children under the age of 5 years worldwide, with long-term consequences on child development and higher risks for several chronic diseases ([@bib46]). Stunting is a multifactorial disease linked to nutritional deficiency, repeated diarrheal infections, poor sanitation, and altered gut bacterial communities ([@bib4], [@bib68]). As phages are key to the maintenance and function of many ecosystems via their regulation of bacterial hosts, we hypothesized they could play a critical role in controlling gut bacterial communities in stunting ([@bib18], [@bib60]). Through a combination of microscopy, 16S rRNA gene sequencing, and metagenomics, we characterized the fecal viral and bacterial communities of Bangladeshi children aged between 14 and 38 months to determine if phage communities differed between non-stunted and stunted children (or undergoing stunting), and identify how phages modify gut bacterial communities of non-stunted and stunted children at a whole community level. We also examined the relative importance of the children's age, diet, and health conditions on bacteria-phage interactions.

Consistent with previous reports, our data indicate high abundances of Proteobacteria and higher abundances of Actinobacteria than expected for their age among the enrolled children, resembling the gut bacterial composition in children younger than 1 year ([@bib23]). Actinobacteria and Proteobacteria are replaced during child development by Firmicutes and Bacteroidetes following the dietary transition from milk to solid food ([@bib43], [@bib55]). Thus, the higher-than-expected levels of Actinobacteria relative to the children's age could suggest a delay in microbiome development and/or maturity in these children. However, these high levels could also be partially explained by the higher prevalence of formula milk in the older children ([@bib23], [@bib49]). Despite higher levels of Proteobacteria in stunted children relative to their non-stunted counterparts, our values are lower than those previously reported in these settings ([@bib43], [@bib49]). Nevertheless, the higher abundances of *Enterobacteriaceae* we observed, known to induce high and persistent levels of intestinal inflammation, coincide with overall lower abundance of Bacteroidetes, widely recognized for their ability to digest complex dietary compounds and improve energy harvested from food ([@bib8], [@bib38]). The higher abundance of *K. pneumoniae* in older stunted children may suggest a role for this taxon in child stunting in this age group. Furthermore, the differences between these age groups were also observed through the distinct bacterial functional traits. Thus, our data are consistent with a developmental immaturity of the gut microbiota in Bangladeshi children.

Using microscopy, we provide the first estimates of free viral particles in stunted and non-stunted children, which are similar for all children irrespective of age and health. These estimates led to VBRs ranging between 0.01 and 0.39, which are lower than those previously reported for nutrient-rich ecosystems ([@bib51]). Although our values fit within estimates for the human gut based on metagenomics ([@bib48]), these low VBRs would support lysogenic replication in the gut. Consistent with findings from previous studies of disease-associated gut phage communities, we observed more reads assigned to the Caudovirales, especially the Siphoviridae family, in stunted children relative to non-stunted ones ([@bib48], [@bib52]). Interestingly, *Lactococcus* phages, which have been associated with Parkinson's disease, were also highly abundant in stunted children younger than 23 months. In contrast, reads assigned to Microviridae and crAssphage were more abundant in non-stunted children. Although both have been identified as members of the gut microbiota of healthy individuals ([@bib34], [@bib61]), their specific roles in human health remain to be determined.

In general, phages followed the same diversity patterns as their bacterial hosts, consistent with the heterogeneity correlation hypothesis ([@bib25], [@bib44]). Indeed, regardless of age, non-stunted children contained high numbers of Firmicutes and Bacteroides phages in their guts, whereas stunted children had lower abundances of Bacteroides phages and higher abundances of Proteobacteria phages (especially in the older children). These convergent diversity patterns were not observed for Actinobacteria and their phages. Armed with this diversity knowledge, we then sought to determine the prevalence of the lytic and lysogenic phage replication cycles, as there is evidence for a switch from lytic to lysogenic replication in the gut with age ([@bib35], [@bib59]). Although the underlying mechanisms remain unknown, it has been suggested that the intrinsic variability in abundance and diversity of gut bacterial communities during child development would support more lytic phages and "kill the winner" dynamics; the more stable gut bacterial communities found in adults, in terms of abundance and diversity, would support lysogenic interactions following "piggyback the winner" dynamics ([@bib35]). However, we report more temperate phages in non-stunted children in both age groups, and highest in the younger children. This suggests that stunted children harbor a gut microbiota and/or specific bacterial metabolism or physiology allowing for lytic phages and kill-the-winner dynamics. This is consistent with previous studies that report a decrease in temperate phages in animal models of IBD and in patients with IBD, where gut bacterial communities are significantly disrupted and altered ([@bib16]). In our case, whether these patterns are due to the local inflamed gut environment, the recurrent diarrheal infections, induction of lysogens, or to specific bacterial taxa that select for lytic phages remains to be determined.

We next proceeded with anaerobic *in vitro* cross-infections of bacterial and viral communities to investigate the role of phage infection in stunting. Our MFA shows a clear clustering of cross-infection conditions distinct from our controls overall and within each age group, demonstrating that these phages are still infectious and actively shaping bacterial community structure *in vitro*. We observed limited changes in bacterial community structure in our heat-killed controls. This could be due to antibacterial effects of phage components and proteins such as endolysins ([@bib54]), the replication of heat-resistant phages, or our incubation conditions that could have triggered prophage induction ([@bib1], [@bib42]). This last observation is supported by our MFA showing the association of heat-killed controls with temperate phages. Alternatively, these changes could be due to labile compounds or peptides, the activity of temperature-resistant small molecules, metabolites, or bacterial debris, which could all selectively promote the growth of a subset of bacterial taxa.

Through cross-infection experiments, we observed age-specific interactions between bacteria and phages communities: phages from children older than 23 months failed to regulate bacterial growth, as seen by increased bacterial abundances. This is further supported by limited changes in bacterial community structure after phage exposure in this age group and differences in major functional traits that contributed the most to β-diversity. In contrast, phages from the younger age group more strongly regulated gut bacterial communities. Specifically, phages from stunted children and children undergoing stunting allowed for Proteobacteria to grow, whereas phages from non-stunted children caused an increase in Firmicutes abundance. Within these younger children, the effects of phage infection on bacterial communities were strongest when bacterial communities from non-stunted children were exposed to phages from stunted children. The changes in phage community composition were due to combined changes in species abundance and turnover. In this age group, *Lactococcus* phages contributed most to the changes in bacterial community diversity of non-stunted children: they were highly abundant in the initial fecal samples of stunted children and strongly correlated with non-stunted bacteria according to our MFA. Their contribution could result from the infection of *Lactococcus* bacteria in non-stunted children, allowing for Proteobacteria (and Enterobacteria) to grow, as suggested by the strong correlation between *Lactococcus* and Enterobacteria in the MFA. There could also be a sex-specific effect of phage communities on bacterial ones, as *Lactococcus* was more prevalent in male children, who also tended to be more breastfed. Sex and breast milk have been associated with distinct gut bacterial and phage communities ([@bib17], [@bib21], [@bib50]), but more data from children with different feeding practices are warranted to determine if this could also affect bacteria-phage interactions.

Bacteria-phage interactions are complex, dynamic, and can be altered by many factors. Recently, it was shown that bacterial and phage abundances can both be maintained at high densities ([@bib9], [@bib10], [@bib59]). The reason why regulation of gut bacterial communities by phages is age-specific remains to be determined. One possible explanation could be the inability of phages to replicate significantly in older children because either they fail to establish a spreading infection or they decline to extinction after briefly propagating. Several mechanisms can affect phage replication, including host availability, defense mechanisms, and phage-phage interactions ([@bib9]). The changes in abundance and diversity of certain bacterial taxa with age could lead to insufficient numbers of adequate hosts to sustain high levels of phage replication ([@bib56]), thereby lowering the efficacy of the phage treatment in the older age group. Alternatively, the differences between these age groups could extend to distinct bacterial defense mechanisms and phage infection properties. It is possible that gut bacterial communities of children older than 23 months, benefiting from a longer exposure to their phage communities and environmental phages, have more diversified bacterial resistance mechanisms resulting from a bacteria-phage arms race coevolution ([@bib25], [@bib59]). This more diversified range of resistance alleles in the bacterial population would limit phage infectivity ([@bib45]), or allow for the emergence of a resistant bacterial subpopulation within a susceptible wild type ([@bib9], [@bib25]). It is also possible that the phages in the gut of younger children are more infectious than in older children and can more effectively target their bacterial hosts ([@bib9], [@bib42]). One must also consider that our experimental setup, using bacterial and phage communities of children from different households, could limit the effect of phages. Indeed, selecting for gut phage communities from non-stunted family members, with higher similarities between bacterial and phage populations, could have increased the efficacy and outcomes of our cross-infections.

Finally, we sought to determine the relative contribution of the children's metadata (diet, age, sex, health, and stunting metrics \[HAZ, WAZ, and WHZ\]) on the changes in bacteria-phage interactions during the cross-infections. Diet, sex, HAZ, WAZ, and WHZ all contributed strongly to shaping bacteria-phage interactions. The effects of these factors on gut bacterial communities have been reported numerous times, and their effect on phage communities might result from the changes in bacterial communities ([@bib21], [@bib66], [@bib68]). Yet, the significant contribution of the Gokushovirinae to bacteria-phage interactions could suggest a direct effect of these variables on phage communities, which remains to be further tested. We found the effects of these factors to be more important for children older than 23 months, most probably due to the weaning process and diversity of solid diets. Our analyses also show distinct differences between the relative contributions of variables associated with younger children from those associated with older ones.

With this study, we determined on a community level that gut phages are distinct in non-stunted and stunted children, and that phage communities isolated from unrelated children can alter the abundance and diversity of gut bacterial communities in an age-specific manner *in vitro*. Importantly, phages from stunted children allowed for Proteobacteria to grow in gut bacterial communities from non-stunted children, hinting at a possible role for phages in the pathophysiology of child stunting. Our pooled sequencing approach precludes the identification of specific bacteria-phage interactions on an individual level but shows community-level interactions that should be further pursued. We observed more active bacteria-phage interactions in early child development, suggesting an intervention time window for microbiome manipulation under the age of 23 months. Going forward, it will be necessary to validate our results *in vivo*, as our data show a strong contribution of host associated factors. It is thus highly likely that other aspects of human health, such as the host immune response and gut inflammation, will also affect bacteria-phage interactions. Using mouse models of stunting ([@bib6]), one could determine if gut phage communities from stunted children enhance the stunting phenotype or pro-inflammatory gut conditions. It will also be crucial to better assess age-specific effects, to refine a possible intervention window, as well as the specific bacteria-phage interactions involved.

STAR★Methods {#sec4}
============

Key Resources Table {#sec4.1}
-------------------

REAGENT or RESOURCESOURCEIDENTIFIER**Biological Samples**Bangladeshi infant stool samplesThis studyN/A**Chemicals, Peptides, and Recombinant Proteins**0.02 um filter, AnodiscWhatman, GE HealthcareCat No./6809-60020.22 um syringe filterMerck-Millipore, USACat No./SLGP033RK0.4 um filter, CycloporeWhatman, GE HealthcareCat No./7060-2504AnaeroGen Gas pack, Oxoid AN0025AThermo ScientificCat No./AN0025ADNase IInvitrogen, Life TechnologiesCat No./AM2222Lysozyme (100 mg/mL)Thermo ScientificCat No./90082Oxygen indicator dry stripBDCat No./271055Protease K (20 mg/ml)Thermo ScientificCat No./AM2548SybrGoldInvitrogen, Life TechnologiesCat No./S11494**Critical Commercial Assays**DNeasy Powersoil KitQiagenCat No./ID: 12888-100Genomic DNA Clean & ConcentratorZYMO researchCat No./ID: D4064illustra Ready-To-Go Genomiphi KitGE healthcareCat No./ID: 25-6601-96MiSeq PE250 instrumentGenome Quebec<http://gqinnovationcenter.com/index.aspx>**Deposited Data**Raw and analyzed data for bacterial and phage metagenomesThis paperSRA: PRJNA594824Raw and analyzed data for bacterial 16S (V4 region)This paperSRA: PRJNA595506**Oligonucleotides**16S rDNA barcodes and primers[Table S3](#mmc4){ref-type="supplementary-material"}[Table S3](#mmc4){ref-type="supplementary-material"}**Software and Algorithms**Bowtie2[@bib27]<http://bowtie-bio.sourceforge.net/bowtie2/index.shtml>BBMergeN/A<https://jgi.doe.gov/data-and-tools/bbtools/>DIAMOND v.0.9.21.122[@bib7]<https://ab.inf.uni-tuebingen.de/software/diamond/>Geneious 11N/A<https://www.geneious.com/>Glimmer[@bib14]<https://ccb.jhu.edu/software/glimmer/>Mauve[@bib13]<http://darlinglab.org/mauve/mauve.html>MIRA[@bib12]<https://sourceforge.net/projects/mira-assembler/>Python programming languagePython Software Foundation<https://www.python.org>Python library MatplotlibN/A<https://matplotlib.org>Python library Numpy v1.12.1N/A<http://www.numpy.org/>Python library Pandas v.0.20.1N/A<https://pandas.pydata.org/>R Statistical Computing SoftwareThe R Foundation<https://www.r-project.org/>adespatial R packageN/A<https://cran.r-project.org/web/packages/adespatial/index.html>FactoMine R packageN/A<http://factominer.free.fr>stats R packageN/A<https://stat.ethz.ch/R-manual/R-devel/library/stats/html/00Index.html>vegan R packageN/A<https://cran.r-project.org/web/packages/vegan/index.html>QIIME2[@bib11]<https://docs.qiime2.org/2019.1/>UCHIME[@bib19]<https://drive5.com/usearch/manual/uchime_algo.html>

Lead Contact and Material Availability {#sec4.2}
--------------------------------------

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Corinne Maurice (<corinne.maurice@mcgill.ca>). This study did not generate new unique reagents.

Experimental Model and Subject Details {#sec4.3}
--------------------------------------

### Study Site and Volunteer Recruitment {#sec4.3.1}

A slum area in Mirpur (location: 23° 48′ 33.51″ N, 90° 21′ 39.34″ E) located in the capital city of Dhaka, Bangladesh was selected as the study site, as detailed elsewhere ([@bib2]). Because diet is a major regulator of the gut microbiota and is tightly linked to child growth, we explored the bacteriaphage dynamics in two groups of children differing in their age-specific diet and having either stunted or normal growth. Anthropometric measurements (height, weight) were taken following a cross-sectional prescreening of children aged between 6-23 months and 24-48 months. These age groups were preliminary targeted because of the strong dietary differences between them. However, upon observation of feeding behavior and recruitment of children at the study site during cross-sectional screening, children enrolled in the younger age group were between 14-23 months old, age-sex matched, and fed a mixture of milk and solid food. Children in the older age group had entirely transitioned to solid food for more than 2 months, were between 28-38 months old, and had a boy-girl ratio of 2:3. Sampling children of both sexes also allowed us to determine if sex was associated with any of our observations. Each group consisted of 15 stunted or undergoing stunting children, and 15 non-stunted children. Nutritional status related to height-for-age Z (HAZ), weight-for-age Z (WAZ) and weight-for-height Z (WHZ) scores were assessed using the [WHO growth reference]{.ul} for the same age and sex (WHO Anthro software, version 3.2.2, WHO, Geneva, Switzerland). The WHO-recommended HAZ cut-off point of \<-2 SD for low height-for-age status was used to define stunting (i.e., when HAZ is below -2 SD of the reference value) in this study ([@bib3], [@bib58]). Children selected for this study had not taken any nutritional supplements, did not receive any antibiotics treatment, nor had a diarrheal episode at least one month prior to sampling. Fecal samples from a total of 60 children (30 stunted and 30 non-stunted) were collected for the study in a pair-wise manner (equal number of stunted and non-stunted samples at a time per group). Sampling was performed between November 2016 and August 2017. This study was conducted following the IRB approved studies: A04M2715B (McGill University) and PR-16001 (icddr,b). Because of the young age of the participants, consent was obtained from their parents/legal guardian, who were provided with oral and written information about the project and research in English and in Bengali. Participation was entirely voluntary.

Method Details {#sec4.4}
--------------

### Sample Collection and Storage {#sec4.4.1}

Immediately after sample collection (\<5 min), a portion of the fresh stool sample was placed inside a 50mL Falcon tube with a gas exchange lid (Bio-reaction tube, Ultident cat\# 229475), which was then rapidly placed inside an anoxic Ziplock bag containing an AnaeroGen Gas pack (Oxoid AN0025A), typically used for the storage of clinical oxygen-sensitive samples. To limit oxygen exposure, a new gas pack was placed inside and the bag was sealed (\<1 min). The absence of oxygen was visually confirmed with an oxygen indicator dry strip (BD BBL). Within 3 hours of sample collection, samples were brought back to the icddr,b, where they were immediately placed in a portable glove box with nitrogen and mixed gas (5% H~2~, 20% CO~2~, balance nitrogen). All subsequent sample-processing steps were performed anaerobically in pre-reduced buffers and media containing cysteine.

### Isolation of Phages and Bacteria from Stool {#sec4.4.2}

Samples from each donor, non-stunted and stunted, were mixed 1 to 10 W:V in reduced Phosphate-buffered saline (rPBS), and vortexed thoroughly to separate phages and bacterial communities from large organic particles matter in the samples. The achieved suspensions were subsequently centrifuged at 700 g for 1 minute, to remove larger organic particles, followed by a 30-minute centrifugation at 3,200g to pellet the bacterial community and separate them from phages. Phage supernatants were transferred to a new tube and bacterial pellets were resuspended in reduced 10% BHI. The phage supernatants were passed through sterile syringe filters (0.22 um, Millex-GP, Merck-Millipore, USA) to separate phages from bacteria that remained in the supernatant. One aliquot of the phage samples from both health status was heat-killed by incubation at 95°C for 20min and 10 min at -80°C followed by DNase (DNase I, Invitrogen, Life Technologies) treatment for 2 h, as a control for the effect of the live phages in the experiment.

### Bacterial and Phage Abundances {#sec4.4.3}

Samples from the same age group were pooled randomly for phage and bacterial enumeration via epifluorescence microscopy (pool size was 3-5 samples). Briefly, phage and bacterial samples were fixed by 1% formaldehyde and diluted in Tris-EDTA buffer. Diluted samples were filtered on 0.4 um (Cyclopore®, Whatman, GE Healthcare) or 0.02 um (Anodisc™, Whatman, GE Healthcare) filters and stained with SybrGold (2.5X final concentration, Invitrogen, Life Technologies) for bacterial and phage enumeration, respectively. Bacterial and phage filters were made in triplicate, except for samples from the younger children, which were done in five replicates. An optical fluorescence microscope (Olympus BX53, Japan) was used to analyze the phage and bacterial filters, and 20 fields of view on average were counted. Bacterial and phage counts were performed on the same day to ensure correct sample dilution for the *in vitro* cross infection experiments.

### *In Vitro* Cross Infections {#sec4.4.4}

A total of 15 experiments were performed on each age group (younger age group: 14-23 months, older age group: 28-38 months). Each of the *in vitro* cross infection experiments was individually performed under anaerobic conditions. On every sampling day, phages and bacterial communities were isolated from fresh fecal samples of one non-stunted and one stunted child via consecutive centrifugation and filtration. The isolated bacterial communities were diluted (∼100x in rPBS) to reach the concentration of ∼10e8.ml^-1^. A similar dilution factor was applied to the phages to ensure bacterial and phage communities were mixed in abundances similar to those determined in the initial stool samples. We proceeded with a full factorial design, where the gut bacterial community from a non-stunted donor was cross-infected with phages from a stunted donor, and vice versa, for 12 hours at 37°C in a diluted nutrient-rich reduced media (10% BHI, brain heart infusion broth) to limit the effect of *in vitro* growth on bacteria-phage interactions. Controls for each cross-infection experiment consisted of (1) non-stunted bacteria with non-stunted phages, (2) stunted bacteria with stunted phages, as well as (3) non-stunted bacteria with heat-killed stunted phages, and (4) stunted bacteria with heat-killed non-stunted phages. Heat-killing effectively reduced the number of phages (free VLPs), which were undetectable at the start of the cross-infection experiments ([Figures 4](#fig4){ref-type="fig"}A and 4B).

### DNA Extraction and Sequencing {#sec4.4.5}

To acquire community-level information and obtain enough material for bacterial and phage metagenomics sequencing, fecal samples were pooled as follows: all the original samples from each age group and for each health status were pooled to generate four pools of samples for phage metagenomics and four pools for bacterial metagenomics. For 16S rDNA sequencing, we created 3 pools of 5 samples per age group per health status, leading to a total of 12 pools. The same strategy was applied to the samples after the *in vitro* cross-infections. We did not find any evidence of the effect of pooling on the bacterial community (phyla and family level) when testing for differences among our three pools with the cross-infection data (family and phyla bacteria), and controlling for sex in the Mixed Regression Tree.

*[Phages.]{.ul}* Fecal samples were first mixed with PBS 1:10 w/v and thoroughly vortexed to separate phage and bacterial communities from organic biological matter. This fecal slurry was then centrifuged at 700xg for 1 minute to remove the larger particles, followed by 30 minutes of 6,000xg centrifugation to separate phages and bacteria. After cross-infections, samples were pooled across treatments and across controls to obtain sufficient genomic material for sequencing. The supernatants were subsequently passed through 0.22 um syringe filters to remove remaining bacteria. The filtrates were centrifuged at 35,000 g for 2 hours to collect the phage particles in the pellets. The phage pellets were resuspended in 1 ml of SM buffer (NaCl/MgSO4⋅7H2O/Tris-Cl/H2O). The supernatants were treated with 40 μL of lysozyme (50 mg.mL^-1^) for 30 min at 37° C; followed by a DNase I treatment for 2h at 37°C to remove non-phage derived DNA. DNase I was inactivated by heating the suspension at 65° C for 30 minutes, and left at room temperature to cool before the next enzymatic reaction. Ten ul of 20% SDS and 40 ul of protease K (20 mg.mL^-1^) were added to each 500ul of sample, and incubated for 1h at 37° C. Subsequently, 35 ul of 5M NaCl and 28 ul of 10% CTAB/0.7 M NaCl were added to the mixture and incubated at 65°C for 30 min. The aqueous phase was transferred to a new tube and 100% cold ethanol was added up to two times of the total volume and incubated overnight at -80°C. Tubes were then centrifuged at 16,000g at 4°C for 1h. Supernatants were removed, and DNA pellets were resuspended in 100 ul of Tris-EDTA buffer. The extracted DNA was purified using the Genomic DNA Clean & Concentrator, ZYMO research (Cat No./ID: D4064), according to the manufacturer's protocol. The concentration of the purified DNA was measured via Qubit, and amplified in triplicate using an Illustra Ready-To-Go Genomiphi Kit, GE, (Cat No./ID: 25-6601-96) according to the manufacturer\'s protocol. Amplification replicates were combined and the isolated phage DNA was sequenced using an Illumina MiSeq PE250 at the McGill University and Genome Quebec Innovation Center (<http://gqinnovationcenter.com/index.aspx>).

*[Bacteria]{.ul}*. Isolated bacterial samples were pooled to make three biological replicates (three replicates per age group and health status). Bacterial DNA was extracted using DNeasy Powersoil Kit, Qiagen (Cat No./ID: 12888-100) according to the manufacturer\'s protocol. The v4 region of the 16S rDNA gene was amplified using triplicate PCR reactions with custom barcoded primers and the pooled amplicons were sequenced as previously detailed ([@bib36]). We also proceeded with shotgun sequencing using an Illumina MiSeq PE250.

### Metagenome Assembly and Analysis {#sec4.4.6}

An overview of the metagenomics analysis is provided [Figure S1](#mmc1){ref-type="supplementary-material"}B. On average 638,978 reads per sample were generated by sequencing for phages; and 789,495 reads per sample for bacteria. For phage analysis, we excluded reads showing nucleotide identity to human genomes (Whole Genome Shotgun, NCBI), representing 0.5% of all reads. A small proportion (2.6%) of the remaining contigs showed nucleotide identity to the NCBI Bacterial databases. Phages are important vectors for horizontal gene transfer, and bacteria often carry one or more prophage as part of their genomes. As such, the removal of all contigs with similarity to bacterial genomes could lead to the loss of phage reads ([@bib41], [@bib52]). We therefore determined whether this small proportion of contigs also showed similarity to phage proteins. Most (86%) of these contigs showed similarity to phage proteins, and were thus included in our analyses. Raw reads produced by sequencing were trimmed for the adapters and quality checked using BBDuk (<http://jgi.doe.gov/data-and-tools/bb-tools/>); default options were used. Trimmed reads were subsequently merged by FLASH([@bib32]). Dedupe (<https://docs.dedupe.io/en/latest/>) from BBTools was used to remove duplicate reads. We used a two-step assembly approach (step 1: De Novo assembly, step 2; Map to reference assembly) for building a cumulative contigs library. MIRA (<https://sourceforge.net/projects/mira-assembler/>) ([@bib12]) was used for the De Novo assembly and Bowtie 2 ([@bib27]) and Geneious 11 (<https://www.geneious.com/>) were used for map to reference assembly (six steps, 3&3)([@bib63]). UCHIME was used to remove chimera ([@bib19]). ORFs were identified for contigs using Glimmer (<https://ccb.jhu.edu/software/glimmer/>) ([@bib14]). All contigs longer than 1,000 bp with a minimum of one ORF were compared to a viral genome reference (downloaded on October 3, 2017, from NCBI and PHAST prophage database) using BlastX (E-value \< 1e-5) for viruses, whereas Megablast was used for bacteria. When reads related to sequences of multiple taxa, a common ancestor was identified. The contigs produced for each health condition were compared using Mauve ([@bib13]) and merged to make a united library to analyze the relative richness of each phage genome. To measure the relative abundance, reads were recruited to contigs generated using Geneious 11 with over 70% coverage and a minimum of 90% identity ([@bib34], [@bib53]). The analysis was performed within Geneious 11. Temperate phages were differentiated from non-temperate based on their nucleotide identity (\>70% coverage and \>90% identity) to known temperate phages (assigned using BlastX, E-value \<1e-5) and/or the presence of lysogeny-associated genes, such as the *integrase* gene. Reads showing \>90% nucleotide identity to an integrase gene reference database (downloaded on March 7, 2019 from NCBI) with over 70% reads' coverage were classified as temperate.

*[E. coli.]{.ul}* While many *E. coli* strains are naturally present in the gut, several pathogenic strains have been characterized, such as Shiga toxin-producers. Contigs that showed nucleotide identity (using a cut-off value of 10e-5, Megablast) to NCBI *E. coli* sequences were determined. The quality-controlled reads were recruited (with \>70% coverage and \>90% identity) to the contigs assigned to different *E. coli* strains.

### Functional Profiles {#sec4.4.7}

An overview of the functional analysis is provided in [Figure S2](#mmc1){ref-type="supplementary-material"}C. Reads were trimmed for adapters and quality checked using BBDuk (<https://jgi.doe.gov/data-and-tools/bbtools/>); a cutoff value of 25 Q score was used. Reads shorter than 60bp were removed. Trimmed reads were merged by BBMerge (<https://jgi.doe.gov/data-and-tools/bbtools/>). Reads subsequently were mapped to Integrated Gene Catalogs (IGC), an integrated catalog of reference genes in the human gut microbiome ([@bib28]), by BLASTX using DIAMOND v.0.9.21.122 ([@bib7])with maximum e-value cutoff 10e3. The best hits provided by IGC were used ([@bib44]).

### 16S rRNA Gene Amplicon Identification {#sec4.4.8}

The V4 hypervariable region was amplified using barcoded 515/806 primer pairs. For the group of younger children, an average of 242,273.7 reads per sample was sequenced (range: 57,964 to 471,403 reads per sample). For the group of older children, an average of 334,750.75 reads per sample was sequenced (range: 269,150 to 420,998 reads per sample). Trimming, merging of paired-end reads, and quality filtering were performed using DADA2. After quality filtering, the group of younger children had an average of 36,789.2 reads per sample, ranging from 11,725 to 63,290 reads per sample with 1683 unique features. For the group of older children after quality filtering, there was an average of 52,775.8 reads per sample, ranging from 35,805 to 72,291 reads per sample with 4,106 unique features. The taxonomic analysis was performed using a pre-trained Naïve Bayes classifier trained on the Greengenes database 13_8 99% OTUs. The 16S rRNA analysis preformed using QIIME2 (<https://docs.qiime2.org/2019.1/>) ([@bib11]).

Quantification and Statistical Analysis {#sec4.5}
---------------------------------------

The comparison of bacterial and virus-like particles absolute abundances between age groups and health status was carried out using Prism 7 (GraphPad software) using a one way test of variance with Holm-Sidak's multiple comparisons post-hoc test. The remaining statistical analyses were conducted in the R open source software (R core team 2018; 36). In particular, the *stats* package was used to run all ANOVAs and t-tests (using the lm() and t.test() functions), the *vegan* package was used to calculate the Shannon diversity (using the diversity() function), *adespatial* was used to obtain the species and treatment contributions to β-diversity (SCBD and LCBD, respectively; using the beta.div() function), and lastly, the *FactoMineR* package was used to run the multiple factor analysis (MFA) (using the function MFA()). Briefly, an MFA is a Principal Component Analysis (PCA) applied to the entire set of variables, but where each set is weighted by their variance in the global analysis. This multivariate correlation analysis is thus useful to examine the relationships among several different types of data or factors, and as such, allowed us to assess the contribution of each factor to changes in the gut microbiome. To restrict the number of variables used in the MFA and help its interpretation, we initially only used the bacterial species, phages species, and functional pathways which dominated the community. This was determined by conducting individuals PCAs on each matrix, and identifying the species or trait that contributed the most to the first two axes of the PCA (using the cleanplot.pca() function in R). For both the β-diversity and MFA on bacteria and phage species data, the species matrices were Hellinger transformed (square root of the Chord transformation) to control for double-zeros (that is, to avoid two observations with 0 abundance to be interpreted as strongly similar in their absence of species). Statistical details of experiments can be found in the Results section and in the figure legends.

Data and Code Availability {#sec4.6}
--------------------------

Metagenomes and bacterial 16S rRNA data can be accessed on the SRA database, accession numbers PRJNA: 594824 (phage and bacterial metagenomes) and PRJNA: 595506 (16S).

Supplemental Information {#app2}
========================

Document S1. Figures S1--S8Table S1. Anthropometric Information of All Children Recruited for this Study, Related to All FiguresSummary of the metadata collected from all 60 children recruited. Categories were transformed as follows: milk, 1 = breastfed, 2 = breastfed + formula milk, 3 = formulae milk; feed, 1 = milk + solid; 2 = only solid; sex, 1 = male, 2 = female; age, month at the day of stool collection.Table S2. Identification of *E. coli* Pathogenic and Non-pathogenic Strains, Related to Figure 2All *E. coli* strains detected from our shotgun data were identified with the following available NCBI references.Table S3. 16S rDNA Barcoded Primer Information, Related to Figure 1Sequencing data (barcode sequences, forward and reverse primer sequences, tag sequences) and specific experimental data (child age group, pool number, sample type, treatment) for the 16S rDNA libraries prepared for this study.Document S2. Article plus Supplemental Information

This work was funded by the Bill and Melinda Gates Foundation (OPP1139814) and the Canada Research Chair Program (950-230748 X-242502) to C.F.M. J.R. and L.D. are funded by the EU ERC StG (803077/PHARMS). The authors thank members of the Maurice lab, J. Shapiro, E. Haggård-ljungquist, H.J. Haiser, and H. Foroughi-Asl for their constructive comments. The graphical abstract was made using Smart Medical Art (SMART): <https://smart.servier.com>.

Author Contributions {#sec5}
====================

M.K.M., D.M., and C.F.M. designed the research; M.K.M., M.A.A.K., P.G., and R.C. performed the experiments; M.K.M., M.T., Z.E.T., M.A.A.K., J.R., L.D., and C.F.M. analyzed the data; M.M.K. oversaw volunteer recruitment and sample coordination; M.K.M., Z.E.T., and C.F.M. wrote the paper. All authors approved this version of the manuscript.

Declaration of Interests {#sec6}
========================

The authors declare no competing interests.

Supplemental Information can be found online at <https://doi.org/10.1016/j.chom.2020.01.004>.

[^1]: Lead Contact
