Ethics
The protocol (registered on ClinicalTrials.gov, study ID number: NCT05629624) was co-developed by researchers from The University of The West Indies, Boston Children’s Hospital, USA, International Centre for Diarrhoeal Disease Research, Bangladesh (icddr,b), and University of Auckland, New Zealand. Ethical approvals were obtained from the Research Review Committee (RRC; August 21, 2021) and Ethical Review Committee (ERC) of icddr,b (protocol no: PR-21084; September 21, 2021), University of Auckland, New Zealand (approval AH23922; for analyses of collected biological samples) and University of West Indies (CREC-MN.51, 21/22).
Study Design and Participants
The study was performed on the baseline data from two cohorts of children who were enrolled (between February – December 2022) as part of the M4EFaD intervention, a community based clinical trial within the Mirpur slum, Dhaka, Bangladesh31. The cohort consisted of a control group of 74 well-nourished children at 12 ± 1 months (WLZ score > -1 SD), an intervention group of 156 children with WLZ < -2 and > -3 and/or MUAC < 12.5 and > 11.5 cm and having MAM at 12 ± 1 months. Inclusion criteria included an assessment of MAM, no history of chronic medical conditions, known congenital anomalies and no antibiotic use within the past month.
Recruitment and anthropometric data collection
Enrolment was initiated on February 7, 2022, and will continue until February 2024. Study surveillance workers (SWs) conducted a door-to-door census (approximately 100,000 households) in Mirpur DNCC wards ward 2, 3 and 5 between January and December 2022. Verbal consent was obtained to participate in the census. The census identified 5736 children aged between 11 to 13 months and 2,314 children aged between 34 to 38 months. During the census, if the guardian verbally consented to the study procedure, and the babies met the inclusion and exclusion criteria of the study, the SWs proceeded to measure the MUAC of the child. Mothers of babies who were within the MUAC range were invited to visit the icddr,b study clinic for further assessment and enrolment.
Final screening for eligibility and study consent occurred at the icddr,b Mirpur study clinic. The consenting process was tailored to each mother’s literacy level and involved reviewing the inclusion and exclusion criteria. Comprehension of the study was assessed using scripted points and open-ended questions.
Following consent, the clinical screening team completed a screening form, capturing the date of enrolment, sex, date of birth (DOB), weight (kg), length/ height (cm), head circumference (cm), and MUAC (cm) measurements of the child. The WLZ/WHZ score for each child was calculated using the WHO anthropometric calculator. The child’s age was validated using the EPI vaccination card. Bayley scores, EEG data were collected upon enrolment to evaluate brain and cognitive development.
EEG data collection and analysis
Continuous scalp EEG was recorded using NetStation 4.5.4. and 128-channel Hydrocel Geodesic Sensor Nets modified to remove eye electrodes (Electrical Geodesics, Inc. (EGI), Eugene, OR, USA)(Fig. S2). Data was sampled at 500 Hz. Impedances were kept under 100 KΩ when possible and measured once at the beginning of the session, and again halfway through. Sessions were conducted in a dimly lit room with the participants sitting on the parent’s lap. The participants were separated from the research staff conducting the session by a curtain, but the testing area was not acoustically or electrically shielded. A second research staff member was present in the testing area to help keep the participant engaged. The subsequent (pre-)processing steps were applied to the resting state data where participants watched a 3-minute video that featured toys.
EEG data were pre-processed with MatLab (R2021B) using the Harvard Automated Processing Pipeline for Electroencephalography (HAPPE) Version 355,56. A specified subset of 30 channels was excluded (E1, E8, E14, E17, E21, E25, E32, E38, E43, E44, E48, E49, E56, E63, E68, E73, E81, E88, E94, E99, E107, E113, E114, E119, E120, E121, E125, E126, E127, E128). Data were bandpass filtered (1-100 Hz), and filtered using a 50 Hz cleanline filter for line noise removal. Bad channels were then automatically identified and rejected, and wavelet-enhanced Independent Component Analysis (ICA) and the Multiple Artifact Rejection Algorithm (MARA) were performed to detect and impute artifacts. Resting state data were segmented into 2s epochs; epochs with an amplitude > ± 150 mV were rejected. Segments were also rejected using segment similarity criteria. Data were then re-referenced to the average of all channels.
EEG outputs from HAPPE were then reformatted and processed using the Batch Electroencephalography Automated Processing Platform (BEAPP)57 to extract power spectra for each participant across the following frequency bands: delta (2–4 Hz), theta (4–6 Hz), low alpha (6–9 Hz), high alpha (9–12 Hz), beta (12–30 Hz), and gamma (30–45 Hz) and the following regions of interest (Fig. S2): occipital (E70, E71, E75, E76, E83), temporal (E36, E40, E41, E45, E46, E102, E103, E104, E108, E109), parietal (E52, E53, E59, E60, E85, E86, E91, E92), and frontal (E5, E6, E12, E13, E24, E27, E28, E33, E34, E112, E116, E117, E122, E123, E124). Further, PSD values were normalized by a log10 transform.
Developmental Outcomes
The Bayley Scales of Infant and Toddler Development, Fourth Edition (BSID-IV) cognitive, language, and motor subscales were administered to all participants35. Research assistants were trained to research reliability in the administration and scoring of the Bayley-4. Due to cultural differences between Bangladesh and the United States where the assessment was developed, Bangladeshi researchers modified some assessment stimuli to improve cultural responsiveness and relevancy. Pictures for the item naming series and action naming series of the expressive language and receptive language subscales were adapted to include items that Bangladeshi children are more likely to be familiar with. For example, in an item depicting a child in a one-piece pyjama set requiring children to predict what the child was about to do, the pyjama set was replaced by bedtime clothing familiar to Bangladeshi children that would signify that the child was going to sleep. The Wolke ratings are a measure of development that have been validated in previous studies of nutritional impacts on behavioural development in Bangladesh36,58–60.
Biological sample collection
Stool samples were collected from each infant at their home at the baseline visit. Samples were collected in DNA/RNA Shield Faecal Collection Tubes (Zymo Research, #R1101). Peripheral venous blood samples were collected in EDTA Vacutainers, separated into plasma and RBCs and immediately frozen at -80 C. Batches of blood and stool samples were air-freighted on dry ice from Bangladesh to the Liggins Institute, New Zealand for processing and analysis.
Microbiome DNA extraction, sequencing, and profiling
DNA was extracted from stool samples using the ZymoBIOMICS MagBead DNA/RNA extraction kit (Zymo Research, #R2136) following the standard protocol. Samples (1 mL) were mechanically lysed in bead bashing tubes using the MiniG tissue homogenizer prior to extraction of DNA. 200 µL of the sample was used post-bead bashing for extraction of DNA following the protocol. A volume of 50 L of elute was collected in DNAse/RNAse Free Water. Samples with a DNA concentration < 14.5 ng/L were re-extracted following the ZymoBIOMICS DNA extraction protocol. Samples were sequenced (Illumina NovaSeq 150PE reads) to an average sequencing depth of 20M read-pairs/sample. Raw sequences were processed using BioBakery3 tools61, specifically read quality filtering and human decontamination with KneadData (Version 1), taxonomic profiling with MetaPhlAn3 (Version 3.1, using the mpa_v31_CHOCOPhlAn_201901 database) and functional profiling using presence/absence and abundance of microbial pathways (MetaCyc) with HUMAnN3 (Version 3.6). A minimum threshold of > 10% prevalence for all detected species was applied for differential abundance analysis.
Plasma lipidomics
Plasma samples for lipidomics were thawed on ice and extracted according to a method modified from62. Briefly, 10 L volume was placed in an amber glass autosampler vial and 300 L of a mixture of Type 1 water, butanol, methanol, chloroform and SPLASH Lipidomix in a ratio of 4:15:15:20:1 was added. The mixture was vortexed and sonicated at room temperature before the protein precipitate was removed by centrifugation and an aliquot of supernatant transferred to an amber glass autosampler for negative ionisation LC-MS/MS. A second aliquot of supernatant was diluted 5 times with 75% IPA for positive ionisation LC-MS/MS. A 5 µL volume of each sample was injected onto a Phenomenex Kinetex F5 column (100 mm × 2.1 mm × 2.6 m) and lipids were separated using a ternary gradient of Type 1 water, methanol and isopropanol containing ammonium acetate. Lipids were quantified and identified with a Q-Exactive mass spectrometer (Thermo Fisher Scientific, Germany) equipped with a heated electrospray ionisation HESI source. Data was processed using MS-DIAL v4.92 9263. For full methodological details see the supplementary information.
Statistical Analyses
Python version 3.9.2 was used to perform all analysis64. Due to the unequal sample sizes and non-normally distributed data; non-parametric statistical approaches were used for differential abundance analysis. Relative abundances were adjusted by Centred Log Ratio to account for the compositional nature of the dataset65. Log adjusted fold change significance was measured using (MWU) test using the ‘mannwhitneyu’ function from ‘scipy.stats’ and adjusted for multiple testing using the ‘fdrcorrection’ function from statsmodels.stats.multitest. PCoA ordinations (plotted using ’skbio.stats.ordination.pcoa’ module) were used to visualise the clustering of the Bray-Curtis dissimilarities (calculated using skbio.distance.pdist) between samples from their species and functional composition. To quantify the variance of the gut microbiome explained covariates, PERMANOVA p-values were calculated from those Bray-Curtis Dissimilarities using the ’permanova’ function from the ’skbio.stats.distance’ module. Associations between features were measured with Spearman correlation (calculated using ’spearmanr’ function from ’scipy.stats’ module), where significance was defined as FDR adjusted p-values66 of < 0.05. Associations between categorical data were measured with Fisher’s Exact test (calculated using ’fisher_exact’ from ’scipy.stats’ module), where significance was defined as p-values of < 0.05.
Machine learning
Random Forest classifiers were trained to differentiate MAM from well-nourished controls. Ensemble models were trained on data that was individually scaled and concatenated. Model hyperparameters including the number of trees in the forest, maximum tree depth, and minimum sample numbers needed to split internal nodes were tuned using grid searching. Model performance was measured with AUCROC. SHAP Value (SHapley Additive exPlanations) interpretation was used to interpret the contributions each feature had on the model’s performance using the ‘shap’ python package67. The AUCROC and mean SHAP Values of a total of 100 Random Forest models for each prediction was used to measure the performance of each hyperparameter combination and to identify overfitting.
Network analysis
Gut bacterial species, gut functional profiles, EEG PSD, Bayley, Wolkes scores, and plasma lipids were included in the network analysis. Significant Spearman rho (q < 0.05) between features that had high mean SHAP scores (SHAP > 85% other values) were used as edges for random forest classifier models that distinguish MAM from well-nourished conditions. Networks were visualised using Cysoscape.