The gut microbiome and early-life growth in a population with high prevalence of stunting

Sub-study population characteristics

The fecal microbiota was characterized in 875 fecal samples from 335 children from 1-18 months of age (Fig. S1). A mean (SD) of 2.6 (1.3) samples were analysed per child. The children in the microbiome sub-study largely resembled the population of all live-born infants in the overall SHINE trial cohort (Table S1), however, the microbiome sub-study included a larger number of children who were born to women living with HIV (29.6%) compared to the whole SHINE cohort (15.6%), due to the deliberate over-sampling of HIV-positive mothers and their infants. In addition, the microbiome sub-study included infants with slightly older mothers and longer gestational ages. The majority of infants were born by vaginal delivery (94.5%) in an institution (89.9%) and were exclusively breastfed (91% at 3 months). Prevalence of stunting (LAZ <-2) varied from 18-34% across study time-points.

Metagenome sequencing performance

Overall, 875 unique infant-visit whole metagenome sequencing datasets were used. On average, 12 million ± 4.2 million quality-filtered read pairs were generated per sample. Sixteen negative controls produced a mean of 655 quality-filtered reads (range = 149 to 1,425; SD = 456). The median percent of human reads detected was 0.05% but ranged widely by age group and decreased over time (Fig. S2a). The median percent un-annotatable reads detected in each sample was 58.6% and increased over time (Fig. S2b). Thirty-six samples were subject to repeated extraction and metagenome sequencing to assess technical variation. These samples originated from 4 unique children, each with 3 visit samples, where each visit sample was extracted and sequenced in 3 replicates. Principal coordinates analysis (PCoA) of Bray-Curtis distances and phylum-level relative abundances revealed little variation between replicates (Fig. S2c-d).

Succession of gut microbiome composition in early childhood

After prevalence and relative abundance threshold filtering, 161 annotated bacterial species were identified. Seven Eukaryotic and 4 Archaeal species were detected in a small proportion of samples, but these did not meet the prevalence thresholds (Table S2). Bifidobacterium longum was the predominant species at all time-points up to 12 months of age. Four other Bifidobacteria speces (B. breve, B. bifidum, B. pseudocatenulatum, and B. kashiwanohense), Escherichia coli, Bacteroides fragilis and Veillonella species were consistently amongst the most abundant species at the earlier time points before being outnumbered by Faecalibacterium prausnitzii and Prevotella copri at 12 and 18 months of age. Taxonomic a-diversity metrics and gene richness tended to decline or remain stable over the first 4-6 months of life, during exclusive breast-feeding, but increased as expected with infant age from 6-18 months of age (Fig. S3a-b), with the introduction of complementary feeds. A large proportion of variation in both compositional (R² = 0.198, P<0.001) and functional 1-diversity (R² = 0.378, P<0.001) was explained by age (Fig. 1a-b and S3 c-d).

• Download figure
• Open in new tab

Fig. 1. Compositional and functional maturation of the gut microbiome of 335 infants from rural Zimbabwe from 1-18 months of age.

PCoA of Bray-Curtis distances of species (a) and metagenomic pathways (b) coloured by age category. PERMANOVA model results are also plotted. The top 20 features and model pseudo-R² from XGBoost models predicting age using species (c) or pathways (d) are ranked by scaled feature importance and relative abundance (0-1) plotted by age (e-f) to visualize taxonomic and functional microbiome succession from 1-18 months of age.

We employed extreme gradient boosting machines (XGBoost), a machine learning approach, to train and test a model of compositional and functional microbiome maturation, with child age as an outcome. Children who were born to HIV-negative mothers, who were non-stunted at 18 months (LAZ > −2) and had at least 2 stool samples collected were used as a ‘training set’, which was then used to predict child age in the remaining samples (test set). Using species composition, the microbiome was highly predictive of child age (Model pseudo-R² = 0.77, MAE = 1.4 months). This ‘microbiota age’ score was also strongly correlated with biological age in the subset of children from the test set who were also non-stunted and born to HIV-negative mothers (pseudo-R² = 0.67). The species most strongly predictive of age included Faecalibacterium prausnitzii, Blautia wexlerae, Prevotella copri, Staphylococcus hominis, Dorea formicigenerans, Bifidobacterium longum, Agathobaculum butyriciproducens, Bifidobacterium bifidum, Bacteroides thetaiotaomicron, Streptococcus vestibularis and Veillonella parvula (Fig. 1c). Metagenomic pathways also predicted age with high accuracy (Model pseudo-R² = 0.68; MAE = 1.5; Fig. 1d). The pathways most strongly predictive of age included methanogenesis from acetate (METH-ACETATE-PWY), multiple nucleotide and amino acid metabolic pathways, including L-tryptophan biosynthesis (TRPSYN-PWY), purine ribonucleosides degradation (PWY0-1296), pyrimidine deoxyribonucleotides de novo biosynthesis I (PWY-7184), L-histidine degradation I (HISDEG-PWY), dTDP-L-rhamnose biosynthesis I (DTDPRHAMSYN-PWY), flavin biosynthesis I (RIBOSYN2-PWY) and nitrate reduction I (DENITRIFICATION-PWY). This ‘metagenome age’ was also highly correlated with age in the subset of children from the test set who were also non-stunted and born to HIV-negative mothers (pseudo-R² = 0.67). Using these models, we created a microbiota-for-age Z-score (MAZ) and metagenome-for-age Z-score (MetAZ), which accounted for variance of microbiota ages with respect to chronological ages at each study visit (see Methods). The top 20 features contributing most strongly to age predictions are plotted in Fig. 1e-f.

WASH and IYCF interventions have little influence on the infant gut microbiome

We have previously reported in the SHINE trial that the WASH intervention had no impact on infant growth, whilst the IYCF intervention increased LAZ scores by 0.16, leading to a 23% reduction in stunting by 18 months of age. We tested whether these randomized interventions impacted any metrics of gut microbiome diversity or maturity in each age group. By performing principal coordinate analysis (PCoA) on Bray-Curtis distances of taxonomic data and functional data we found no significant differences in 1-diversity between IYCF and non-IYCF arms at any time-point. There was a significant difference in Bray-Curtis distances for microbiome composition between WASH and non-WASH arms at the 3-month time-point (PERMANOVA, P=0.013; R²= 0.011), but at no other time-points (Fig. 2a-d). No significant differences were observed in a-diversity metrics or gene richness between WASH versus non-WASH arms nor IYCF versus non-IYCF arms at any time-point (Fig. 2e-f). Using multivariable regression analyses (MaAsLin2), there were also no differences in the relative abundance of species or pathways between intervention arms apart from a small number of features at 3 months in the WASH arms (increased Klebsiella pneumoniae, reduced Collinsella aerofaciens and more abundant metagenomic pathways involved in biotin and folate synthesis) and at 18 months in the IYCF arms (reduced Eubacterium siraeum, E. rectale and Agathobaculum butyriciproducens; Table S3 and S4).

• Download figure
• Open in new tab

Fig. 2. Impact of randomized WASH and IYCF interventions on infant gut microbiome.

PCoA of Bray-Curtis distances species coloured by WASH vs non-WASH arms (a), including PERMANOVA model results, and IYCF vs non-IYCF arms (b) are plotted in addition to the first component (PC1) from PCoA of species (c) and pathways (d). The IYCF intervention was introduced after 6 months of age, therefore direct comparisons of IYCF vs non-IYCF arms are not shown in the 1-, 2- and 3-month age categories. No significant differences were observed in Shannon alpha diversity (e) and gene richness (f) according to trial arm.

HIV exposure comprehensively alters infant gut microbiome composition and function

We previously reported in this cohort that maternal HIV exposure significantly impacts infant growth, whereby CHEU displayed significantly poorer linear growth compared with children who are HIV-unexposed (CHU) (30). We assessed diversity metrics and microbiome maturity in CHEU versus CHU and found that CHEU displayed significantly greater diversity (Shannon index; Wilcoxon rank-sum test P = 0.002; Fig. 3a) and species richness (P = 0.01; Fig. 3b) compared with CHU at 12 months of age, although the size of the two groups was imbalanced within this age category (n=91 CHEU, n=27 CHU). Other metrics of a-diversity were significantly lower in CHEU at 2 months of age (evenness P = 0.04, Simpson index P = 0.04). Metagenomic gene richness was also elevated in CHEU at 1, 3, 6 and 12 months (Wilcoxon rank-sum test P = 0.002, P = 0.059, P = 0.007, P < 0.001, respectively; Fig. 3c). Analysis of taxonomic 1-diversity between CHEU and CHU also revealed significant differences at 1, 3, 6, 12 and 18 months suggesting that in utero HIV exposure was significantly associated with gut microbiome succession and development throughout the first 18 months of life (Fig. 3d-e). HIV exposure also explained a significant proportion of the variation in metagenome pathway 1-diversity (PERMANOVA, P = 0.002, R2 = 0.033) at 1 month of age (Fig. 3h-i), but not at later ages. We next tested the association between HIV exposure and gut microbiome maturity. We found that CHEU displayed greater microbiome age and MAZ, and hence microbiota over-maturity, compared with CHU at 12 months of age (median 14.7 vs 9.2 months; 1 = 3.18, P < 0.001), which tended to be higher in CHEU at 1 months and at 6 months also (P = 0.074 and P = 0.059, respectively; Fig. 3 f-g). However, at 18 months CHEU displayed lower microbiome age (15.3 vs 16.4 months; 1 = −1.1, P = 0.045 respectively). Similarly, CHEU displayed significantly greater metagenome ages and MetAZ scores compared with CHU at 1 month (1 = 1.1, p = 0.039) and 12 months of age (1 = 2.2, P = 0.04; Fig. 3j-k), suggesting that HIV exposure drives both compositional and functional microbiome over-maturity. The proportion of CHEU receiving prophylactic cotrimoxazole ranged from 56-83% across the 3, 6, 12 and 18-month study visits.

• Download figure
• Open in new tab

Fig. 3. Maternal HIV infection comprehensively alters infant gut microbiome diversity and maturity.

Shannon alpha diversity (a), species richness (b) and gene richness (c) shows significant over-diversification in CHEU vs CHU (Wilcoxon rank-sum test; *p<0.05). PCoA of Bray-Curtis distances (d) and PC1 (e) of species composition in CHEU vs CHU show significant differences throughout 18 months of life (PERMANOVA). Microbiome age (f) and microbiome-for-age Z score (MAZ; g) shows significant differences in gut microbiome maturity in CHEU vs CHU (linear regression analyses). PCoA of microbiome gene pathways shows differences in CHEU vs CHU at 1 month of age (h-i) in addition to differences in metagenomic maturity (j-k).

We next explored which species were differentially abundant between CHEU and CHU by performing multivariable regression analyses, adjusting for age at stool sample collection, exclusive breastfeeding status, delivery mode, and randomised trial arm. Between 1-3 months of age, two Bifidobacteria species, B. longum and B. bifidum (Fig. 4a-b), in addition to Veillonella seminalis were significantly less abundant in CHEU versus CHU (q < 0.1; Table S3). Conversely, Flavonifractor plautii was significantly more abundant in CHEU at 3 months (q = 0.02). At 18 months, B. bifidum was again significantly less abundant in CHEU (q = 0.04), whilst two other Bifidobacteria species were also weakly associated with CHEU, whereby B. breve was lower and B. pseudocatenulatum higher (q < 0.25). Regression analyses of metagenomic pathways with CHEU status generated similar outcomes. Following adjustment for covariates, significant associations for metagenomic pathways were only present at 1 and 3 months of age. At 1 month of age, these included significant negative associations between CHEU and amino acid synthetic pathways (superpathway of L-threonine biosynthesis, superpathway of L-isoleucine biosynthesis I and L-lysine biosynthesis I; Fig. 4c-e) and positive associations with pathways involved in the degradation of sugar derivatives, including fructuronate, glucoronate and galacturonate (PWY 7242 D-fructuronate degradation, PWY 6507 4-deoxy-L-threo-hex-4-enopyranuronate degradation, GALACTUROCAT PWY D-galacturonate degradation I, GALACT GLUCUROCAT PWY superpathway of hexuronide and hexuronate degradation and GLUCUROCAT PWY superpathway of beta D-glucuronide and D-glucuronate degradation; Fig. 4f-h, Table S4). A handful of pathways involved in fatty acid oxidation and fermentation (fatty acid beta oxidation peroxisome and succinate fermentation to butanoate) were also significantly positively associated with CHEU at 3 months of age.

• Download figure
• Open in new tab

Fig. 4. Maternal HIV infection is associated with reduced abundance of Bifidobacteria abundance and amino acid biosynthesis genes.

Relative abundance (0-1) of Bifidobacterium longum (a) and B. bifidum (b) in the gut microbiome of CHEU and CHU at each age category via multivariate regression analyses. Multivariate regression of gene pathways demonstrates reduced abundance of amino acid biosynthetic pathways (c-e) and increase in abundance of pathways involved in degradation of sugar derivatives (f-h). Multivariate regression using MaAsLin2 using default q < 0.25 cut-off. *q<0.1.

Microbiome functionality, but not composition, predicts linear and ponderal growth

We next examined the relationship between taxonomic and functional features of the gut microbiome and attained growth (LAZ and WHZ) and growth velocity (WHZ and LAZ velocity, increase in z-score increments per day between visits) from 1-18 months of age using XGBoost models. Models were run separately for each of the six age groups, stratified by maternal HIV status and run in two combinations of predictive features: (i) microbiome features alone (species or pathways); and (ii) microbiome features and epidemiological variables, which included maternal anthropometry, baseline WASH and infant diet variables, amongst others (Table S5). In models combining microbiome features with epidemiological features, birthweight, maternal height, maternal mid-upper arm circumference, and household wealth were all important predictors of attained infant LAZ and WHZ and growth velocity. Models including microbiome taxonomic features (species) alone performed poorly for both attained and growth velocity at every age category and regardless of HIV exposure status (Fig. 5a and 6a), with a majority of models resulting in pseudo-R² values < 0. Model performance for linear growth improved when epidemiological features were included, suggesting that gut microbiota composition alone was poorly predictive of growth. Taxonomic features were weakly predictive of WHZ velocity at 1-3 months (pseudo-R² 0.09-0.12) and attained WHZ at 6 months (pseudo-R² = 0.22), but only in children born to HIV-negative mothers (Fig. 6a). Conversely, models containing functional metagenomic pathways were moderately predictive of both attained and future growth throughout 18 months of age (pseudo-R² = 0-0.55; Fig. 5a and 6a) albeit with relatively large mean absolute errors (MAE) for both linear (0.48-0.95 LAZ) and ponderal growth models (0.6-1.18 WHZ). The inclusion of epidemiological variables in the metagenomic models added little to performance suggesting that pathway features were independently predictive of both linear and ponderal growth. Models predicted WHZ better than LAZ and models including children born to HIV-negative mothers also tended to perform better. MAE decreased in all models as age increased (Fig. S4 a-d).

• Download figure
• Open in new tab

Fig. 5. Prediction of attained LAZ and LAZ velocity using XGBoost models.

Performance of XGBoost models as assessed by pseudo-R² values for prediction of attained LAZ and LAZ velocity (LAZ increase per day to next study visit) using species or metagenomic pathways, stratified by age category and maternal HIV status (a). Models were run using microbiome features alone (species or metagenomic pathways; blue points) and in combination with epidemiological variables (yellow points). The top ranked pathways predicting LAZ at each age category are plotted (b), stratified by maternal HIV status and coloured by scaled importance in the XGBoost model. Accumulated effect plots (ALE) of representative pathways ranking highly in XGBoost model predictions display change in predicted linear growth (LAZ or LAZ velocity) by percentile of the feature abundance distribution. Tick marks on the x-axis are a rug plot of individual feature abundance percentiles. ALEs were generated using the ALEplot package and were plotted using ggplot2. Standard deviations (sd) were calculated per increment in microbiome feature and were used to calculate and plot increment-wise 95% confidence intervals as the average change in the outcome ±1.96(sd/sqrt(n)), where n is the number of observed feature values, and sd is the standard deviation of the change in the outcome variable in an interval.

• Download figure
• Open in new tab

Fig. 6. Prediction of attained LAZ and LAZ velocity using XGBoost models.

Performance of XGBoost models as assessed by pseudo-R² values for prediction of attained WHZ and WHZ velocity (WHZ increase per day to next study visit) using species or metagenomic pathways, stratified by age category and maternal HIV status (a). Models were run using microbiome features alone (species or metagenomic pathways; blue points) and in combination with epidemiological variables (yellow points). The top ranked pathways predicting WHZ at each age category are plotted (b), stratified by maternal HIV status and coloured by scaled importance in the XGBoost model. Accumulated effect plots (ALE) of representative pathways ranking highly in XGBoost model predictions display change in predicted linear growth (LAZ or LAZ velocity) by percentile of the feature abundance distribution. Tick marks on the x-axis are a rug plot of individual feature abundance percentiles. ALEs were generated using the ALEplot package and were plotted using ggplot2. Standard deviations (sd) were calculated per increment in microbiome feature and were used to calculate and plot increment-wise 95% confidence intervals as the average change in the outcome ±1.96(sd/sqrt(n)), where n is the number of observed feature values, and sd is the standard deviation of the change in the outcome variable in an interval.

Microbiome features associated with linear growth

In all but the 2-month age group, birthweight contributed most strongly to prediction of attained LAZ. Metagenomic pathways consistently performed better as predictors than other epidemiological features including maternal height. The most predictive pathways were largely similar between infants born to HIV-positive and HIV-negative mothers. At 1 month and 2-months, metagenomic pathways encoding enzyme co-factor biosynthesis, nucleotide degradation and biosynthesis and amino acid biosynthesis were consistently predictive of both attained LAZ and LAZ velocity (Fig. 5b). Accumulated local effects (ALE) (36, 37) plots show the average effect of some of the most important features on model outcomes (Fig. 6c). At 3 months and 6 months, pathways encoding fermentation and carbohydrate biosynthesis were consistently predictive of attained LAZ and LAZ velocity, whilst amino acid degradation pathways, amongst others, were predictive of growth at the oldest age groups. In particular, pathways involved in vitamin B biosynthesis (flavin, folate, biotin, thiazole and cobalamin biosynthetic pathways) were consistently predictive of attained LAZ, and included flavin biosynthesis I, 6-hydroxymethyl-dihydropterin diphosphate biosynthesis, superpathway of tetrahydrofolate biosynthesis, adenosylcobalamin salvage from cobinamide I, biotin biosynthesis II and thiazole biosynthesis I. Increasing abundances of B vitamin biosynthesizing genes contributed to increasing predicted growth in ALE plots, apart from Thiamin biosynthesis, which predicted lower LAZ at 12 months (Fig. 5c). At 12 months of age, 4-coumarate degradation (anaerobic), a pathway involved in plant polysaccharide degradation, was the most predictive pathway of LAZ in children born to HIV-positive mothers, with greater abundance associated with greater LAZ. We also assessed pathways predicting growth velocity and found similar results to that of attained growth (Fig. S5a). At 2 months of age, folate biosynthetic pathways (folate transformations II and N10-formyl-tetrahydrofolate biosynthesis) were the top two predictive features of linear growth velocity in children born to HIV-negative mothers, whilst at 3 months, purine and pyrimidine pathways were amongst the most predictive features of linear growth velocity. Similarly, as per the attained growth models, amino acid and fatty acid biosynthetic pathways were also strongly predictive. Glycogen biosynthesis pathways were also consistently predictive of linear growth velocity at all ages from 2 months onwards.

Microbiome features associated with ponderal growth

In the few models incorporating taxonomic features that weakly predicted WHZ velocity in infants born to HIV-negative mothers, Escherichia coli at 2 months and Bacteroides fragilis and Veillonella atypica at 3 months were amongst the most predictive features (Fig. S6). Many of the same categories of biosynthetic microbiota pathways were predictive of both WHZ and LAZ, including amino acid and nucleotide (especially purine biosynthesis) biosynthetic pathways in addition to a number of lipid synthesis pathways at older age groups (Fig. 6b-c). The most predictive pathways were largely similar between infants born to HIV-positive and HIV-negative mothers. O-antigen biosynthesis pathways (PWY-7328 UDP-glucose-derived O-antigen building blocks biosynthesis & PWY-7332 UDP-N-acetylglucosamine-derived O-antigen building blocks biosynthesis and OANTIGEN-PWY pathway), which did not appear in LAZ models, were consistently amongst the most predictive features of WHZ at 1, 6, 12 and 18 months, whereby greater abundance was associated with reduced growth. Pyrimidine and purine synthetic pathways were consistently the strongest predictors of WHZ velocity with varying directions of association, including superpathway of pyrimidine deoxyribonucleotides de novo biosynthesis (PWY0-166) which was the strongest predictive feature of WHZ velocity at 12 months in children born to HIV-mothers (Fig. S5b). Similar to the attained WHZ models, O-antigen biosynthesis pathways, amino acid synthetic pathways and glycogen biosynthesis pathways were all strongly predictive of ponderal growth velocity.

SHINE trial design

The study design and methods for The Sanitation Hygiene Infant Nutrition Efficacy (SHINE) trial and for the corresponding microbiome analyses, have been reported previously (57, 58). Briefly, SHINE was a 2×2 cluster-randomized trial, conducted between 2012 and 2017, to determine the independent and combined effects of improved infant and young child feeding (IYCF) and WASH on child stunting and anaemia in two rural Zimbabwean districts (NCT01824940). 5280 pregnant women were cluster-randomized to one of four interventions: WASH, IYCF, WASH+IYCF, and Standard of Care (SOC). The SOC interventions, included in all trial arms, comprised exclusive breastfeeding promotion for all infants up to 6 months and strengthened prevention of mother to child transmission (PMTCT) of HIV services. The household WASH intervention was initiated during pregnancy and was designed to reduce exposure to human and animal feces, including, at the household level: construction of a ventilated improved pit latrine, installation of two hand-washing stations plus monthly delivery of liquid soap and water chlorination solution, provision of a play space for the infant, and hygiene counseling. The IYCF intervention was designed to improve infant diets using a small-quantity lipid-based nutrient supplement (SQ-LNS), provided to the infant from 6-18 months, and educational interventions promoting the use of age-appropriate, locally available foods and dietary diversity. Lastly, a combined trial arm, WASH+IYCF, evaluated the effects of both improved WASH and infant nutrition.

Infants were followed up at study visits at 1, 3, 6, 12 and 18 months of age. Length and weight were measured at each infant visit, as described previously (35). Length-for-age z scores (LAZ) and weight-for-height z scores (WHZ) were calculated from length and weight measurements at each visit according to WHO Child Growth Standards. Epidemiologic data for the infants was collected from the baseline and follow-up visits using trial questionnaires that included maternal anthropometry, birth outcomes, baseline household WASH facilities, household wealth, maternal education, religion, parity, household size, dietary diversity, changes in breastfeeding and complementary feeding practices, food security, 7-day and 3-month infant health status, and antimicrobial use.

HIV testing

HIV testing was conducted on mothers at the baseline visit using a rapid test algorithm (Alere Determine HIV1/2 test, followed by INSTI HIV-1/2 test if positive). Those testing positive for HIV had CD4 counts measured (Alere Pima Analyser) and referral to local clinics; women were encouraged to begin co-trimoxazole prophylaxis and ART, to exclusively breastfeed, and to attend clinic at 6 weeks postpartum for early infant diagnosis and infant co-trimoxazole prophylaxis. Women testing negative for HIV were offered retesting at 32 gestational weeks and 18 months postpartum. Children of mothers living with HIV were offered testing for HIV at each of the study visits. Those who tested positive were referred to local clinics for ART. HIV was diagnosed using DNA PCR on dried blood-spot samples or RNA PCR on plasma in samples collected prior to 18 months. In samples collected after 18 months, HIV was diagnosed by PCR or rapid test algorithm, depending on samples provided. Children born to HIV-positive mothers and whom tested negative at 18 months were classified as HIV-exposed uninfected (CHEU). Inconclusive or discordant results were re-tested; if no further samples were available or repeat testing was inconclusive, children were classified as HIV-unknown.

Microbiome sub-study

All CHEU and a subgroup of CHU from the SHINE study were enrolled into an Environmental Enteric Dysfunction (EED) sub-study (n=1,656 mother-child pairs); these infants underwent intensive biological specimen collection at 1, 3, 6, 12 and 18 months of age (58). The EED sub-study was therefore enriched for mothers living with HIV, by design. Sample selection for inclusion into the current microbiome study was conducted to enhance longitudinal profiling of the mother and infants gut microbiota. Of the mother-infant pairs within the EED sub-study, those with least one maternal fecal specimen (of 2 possible) and at least 2 infant fecal specimens (of 5 possible) were included in the gut microbiome analyses. An additional 94 samples collected at the 1 and 3-month visits, that did not meet these criteria, but had microbiome sequencing data available from a separate study examining rotavirus vaccine immunogenicity in the SHINE trial (59), were also included in these analyses. Infant ages varied at each study visit due to the allowable window around the visit date for the larger SHINE trial. Therefore, for this microbiome study, each stool samples were re-categorized into 6 age groups corresponding to important stages in infant microbiome development: “1 month” (0-6 weeks), “2 months” (7 weeks – <3 months), “3 months” (3-6 months), “6 months” (6-9 months), “12 months” (9-15 months), and “18 months” (15-20 months).

Sample collection

Study visits were conducted by trained study nurses in participants’ homes. Sterile stool collection tubes were provided to mothers, who collected stool samples from their infants on the morning of each study visit. Samples were placed in cool boxes immediately upon collection by study nurses and transported by motorbike to field laboratories where they were aliquoted and stored at −80°C within 6 hours of collection before subsequent transport to the central laboratory in Harare for long-term storage at −80°C. An aliquot of each stool sample was shipped on dry ice by courier to the British Columbia Centre for Disease Control in Vancouver, Canada. A strict cold chain was maintained throughout transport, ensuring no freeze-thaw cycles occurred between sample collection and processing.

Whole metagenome library preparation and sequencing

DNA was extracted from 100-200mg of stool samples using the Qiagen DNeasy PowerSoil Kit as per the manufacturer’s instructions. DNA quantity was assessed by fluorometry (QuBit) and quality confirmed by spectrophotometry (SimpliNano). 1µg DNA was subsequently used as input for metagenomic sequencing library preparation using the Illumina TruSeq PCR-free library preparation protocol, using custom end-repair, adenylation and ligation enzyme premixes (New England Biolabs). The concentration and size of constructed libraries were assessed by qPCR and by TapeStation (Agilent). DNA-free negative controls and positive controls (ZymoBIOMICS) were included in all DNA extraction and library preparation steps. Libraries were pooled in random batches of 48 samples including one negative control. A set of specimens were subject to replicate DNA extraction, library preparation and, sequencing to estimate the magnitude of technical variability among samples. Whole metagenome sequencing was performed with 125-nucleotide paired-end reads using either the Illumina HiSeq 2500 or HiSeqX platforms at Canada’s Michael Smith Genome Sciences Centre, Vancouver, Canada.

Bioinformatics

Sequenced reads were trimmed of adapters and filtered to remove low-quality, short (<70% raw read length), and duplicate reads, as well as those of human, other animal or plant origin, using KneadData with default settings. Species composition was determined by identifying clade-specific markers from reads using MetaPhlAn3 with default settings (60). Relative abundance estimates were obtained from known assigned reads, and unknown read proportions were estimated from total, assigned and unassigned, reads. Percent human DNA was estimated from KneadData output, using the proportion of quality-filtered reads that align to the human genome. Given the smaller viral genome sizes, sequencing depth, and limitations of MetaPhlAn3 for virus identification, we did not include viruses in our current analyses. We applied a minimum threshold of >0.1% relative abundance and ≥5% prevalence for all detected species. Metabolic pathway composition was determined using HUMAnN3 with default settings against the UniRef90 database (60). Pathway abundance estimates were normalized using reads per kilobase per million mapped reads (RPKM) and then re-normalized to relative abundance. We applied a minimum relative abundance threshold of 3×10^-7% and ≥5% prevalence for all metagenomic pathways.

Statistical analysis

All data were analysed using R (v.4.0.5). Microbiome data were handled using the phyloseq package. Alpha diversity metrics were calculated using the vegan package. Beta-diversity was estimated using the Bray-Curtis dissimilarity index and analysed by permutation analysis of variance (PERMANOVA). Differential abundance analysis of species or functional pathways was assessed using multiple regression analyses using the MaAsLin2 package (61). Four covariates were chosen for adjustment in MaAsLin2 regression models and included age at stool sample collection, exclusive breastfeeding status (recorded at 3 months old), delivery mode, and randomised trial arm. These covariates were chosen based on biological plausibility and previous evidence of their influence on gut microbiome composition in large birth cohorts (17). Adjustment for multiple comparisons was performed using the Benjamini-Hochberg false discovery rate (FDR).

The SHINE trial did not observe an interaction between the randomized WASH and IYCF interventions and growth; therefore, randomised trial arms were combined into WASH versus non-WASH arms and IYCF versus non-IYCF arms for specified analyses. We restricted the IYCF analysis to the 6, 12 and 18 month visits, corresponding to the period during which supplemental infant feeding was introduced (from 6 months of age). All children, regardless of HIV status or exposure status were included in the growth analyses (875 total stool samples). 16 samples had missing ages and were excluded from age prediction models. Stool samples collected from children classified as HIV-unknown (24 samples) or HIV-positive (4 samples) at 18 months were excluded from direct comparisons of CHEU vs CHU infants.

XGBoost models

Relationships between the infant microbiome and age or growth (attained LAZ/WHZ and WHZ/LAZ velocity) were evaluated using extreme gradient boosting machines (XGBoost). XGBoost builds an optimized predictive model by creating an ensemble from a series of weakly predictive models. XGBoost is also non-parametric, can capture non-linear relationships, and can accommodate high-dimensional data (62). The XGBoost models were developed using microbiome relative abundances (species or pathways) and all epidemiologic variables.

XGBoost model selection was performed in 3 stages as previously described (27). In brief, BayesianOptimization function of the rBayesianOptimization package was used with 10-fold cross-validation to select model hyperparameters by minimizing the mean squared error (MSE). Models with the lowest MSE (in the 5^th percentile) were retained, and from these models the variables that contributed to the top 95% of variable importance by proportion were retained. In stage two, all epidemiologic variables were included with the microbiome variables obtained in stage one, and BayesianOptimization was re-run as for stage one but using leave-one-out cross-validation. Microbiome variables that contributed to the top 95% of variable importance by proportion were retained. In stage three, all epidemiologic variables, microbiome features, and hyperparameters selected in stage two were used to fit our final models, using leave-one-out cross-validation to minimize the MSE. This 3-stage hyperparameter tuning and model building was performed for two feature sets, one comprising microbiome features and a second comprising microbiome plus epidemiologic features; this was done to assess model performance and to examine the contribution of epidemiologic versus microbiome features. Separate models were built for attained LAZ and WHZ and growth velocity outcomes (13). We assessed microbiota composition and functional pathways separately. XGBoost models were fit using the H20.ai engine and h2o R package interface with the XGBoost package. XGBoost model performance was evaluated using pseudo-R² and mean absolute error (MAE). Pseudo-R² values < 0 indicated that the prediction of the model was worse than the mean response. Scaled relative importance for each model feature was used to identify the twenty most informative variables for further interpretation, where the most important variable is ranked first, and the importance of subsequent variables are relative to the first variable. The marginal relationships between the twenty most important features and each growth outcome were visualized for interpretation (36) using accumulated local effects plots (ALE). ALE plots can be interpreted as showing a marginal effect, adjusted for all covariates retained in the final model, showing the expected change in the outcome variable per increment in a model feature. The resulting effect sizes are plotted cumulatively and centered about the average effect size (37). ALEs were generated using the ALEplot package, modified to compute confidence intervals, and were plotted using ggplot2. Standard deviations were calculated per increment and were used to calculate and plot increment-wise 95% confidence intervals.

Microbiome age

We investigated microbiome maturation by building an age prediction model using XGBoost and microbiome features only (species or pathways). We partitioned the infants and their corresponding datasets into three groups to train and test a model of microbiome age: (1) CHU with LAZ > −2 at 18 month of age, who contributed >1 dataset (healthy training set), (2) remaining CHU infants with LAZ > −2 at 18 month of age, who contributed a single metagenomic dataset (healthy test set); and (3) CHEU or children with LAZ ≤ −2 at 18 months of age (unhealthy test set). Age was log transformed as a response in the XGBoost model. We performed the same 3-stage tuning and model building procedure, as described above. We generated model performance metrics, including pseudo-R2, mean absolute error (MAE) and mean squared error (MSE) for the three sets. For the training set, we used the cross-validation, hold-out predictions to generate the metrics, and for the two test sets, we used the predicted values from the final to calculate the model performance metrics. We exponentiated the predicted log transformed ages and plotted these values against the observed age. The predicted age using these models is referred to as ‘microbiota age’ for the models trained using species and ‘metagenome age’ for models trained using pathways. To account for variance of microbiota ages with respect to chronological age within the age range of each study visit, a microbiota for age Z-score (MAZ) and metagenome-for-age Z-score (MetAZ) was also created using the microbiome age and metagenome ages as previously described (19). A Z-score was calculated to account for variation in ages within each study visit using the following formula: (Microbiota age of child – median microbiota age of ‘healthy’ child at same study visit)/ standard deviation of microbiota age of ‘healthy’ child at same study visit.