Experimental design

Our a priori hypothesis was that diet impacts the gut microbiota of preterm infants and that infants fed formula would have a higher ARG load than infants fed only human milk, due to the changes in the community composition. To study the effect of formula feeding, we drafted a statistical analysis pipeline, including several steps to ensure model quality and the robustness of the results (Figure 2).

Open in a separate window

FIGURE 2

Flowchart of analysis methods, describing the included cohorts and the statistical analysis pipeline devised for this study. The pipeline was used to ensure that the effect of diet could be validated in several points of the analysis and that suitable distributions and model types are applied to the data. Blue boxes describe steps in preliminary analysis with the training set of 46 infants, and green boxes steps taken in model testing and cross-validation with meta-analysis cohorts. Abbreviations: ARG, antibiotic resistance gene.

As the first step in the statistical analysis pipeline outlined in Figure 2, we used a random forest (RF) machine-learning algorithm to identify explanatory variables linked to ARG load in the preterm infant feces for ARG load model training. We used RF variable importance to prefilter clinical parameters and metadata to select which clinical variables we should use in the subsequent model building.

Secondly, once we determined that the responses were linear, we used gamma-distributed generalized linear models. We chose this distribution because the data were ARG counts normalized with 16S rRNA gene count, and thus were unsuitable for negative binomial or quasi-Poisson distributions, and because the data had overdispersed variances violating the normality assumption. GLMs are a potent tool for hypothesis testing but thus far have been used in only a limited number of microbiome studies (4, 19). After the GLM distribution was selected, we used the variables that could potentially affect ARG load to train the model using the preterm cohort. Explanatory variables that did not significantly explain variance in the ARG loads or did not improve the model were dropped 1 at a time from the full model using Akaike information criterion (AIC) values and chi-squared tests (α = 0.05). Model training ultimately resulted in the simplest model possible explaining the variation in the data without overfitting the model (20).

Thirdly, we confirmed the applicability of the model built using data gathered for a meta-analysis as a test set. Lastly, we created a de novo model, following similar steps as for the train set, to confirm the association of formula with an increased ARG load even when we considered all data and all clinical variables.

Ethical approval

The Institutional Review Board of Pennsylvania State University (IRB #35925) approved the study.

Sample collection

Enrolment commenced in January 2014 and concluded in July 2017. Inclusion criteria were preterm infants born 26–37 weeks of gestation who were admitted to the Penn State Hershey NICU or transferred to the PSHMC NICU within 72 hours of birth. Exclusion criteria included the following: infants born at <26 or >37 weeks of gestation, infants born with major congenital anomalies (heart, gastrointestinal, renal, or respiratory tract), mothers known to use illicit drugs or abuse alcohol, or mothers with a history of depression requiring long-term psychotropic medication. We obtained written consent from all subjects’ parents within 48 hours of admission to the PSHMC NICU. Metadata on feed history, clinical course, and necrotizing enterocolitis (NEC) outcomes were collected electronically during the entire course of stay in the PSHMC NICU.

Fecal material was collected into sterile microcentrifuge tubes approximately 2 weeks after prophylactic antibiotics had been discontinued and enteral feeds were initiated. Samples were frozen at −80°C until analysis.

DNA isolation and quantification

Fecal samples were collected and shipped to Wright Labs, LLC. Nucleic acid extractions were performed from approximately 0.25 g of feces using the Qiagen DNeasy PowerSoil DNA Isolation kit following the manufacturer's instructions (Qiagen). The lysing was performed using the Disruptor Genie cell disruptor (Scientific Industries). Genomic DNA was eluted in 50 μL of 10 mM Tris. Subsequent quantification was performed using the Qubit 2.0 Fluorometer (Life Technologies) using the double-stranded DNA high-sensitivity assay.

Metagenomic sequencing

The Nextera XT library preparation kit (Illumina) with the manufacturer's standard protocol was used for library preparation. The library was paired-end sequenced using 1 NextSeq run targeting an average of 10 million reads per sample. Sequencing and library preparation were done at the Institute of Biotechnology's sequencing services (University of Helsinki).

Metagenomic analysis

We performed quality control using the FastQC and MultiQC programs (21, 22). The sequences were trimmed using Cutadapt version 1.10 with the options -m 1, -e 0.2, -O 15, and -q 20 to filter out adapters and low-quality reads (23). Filtered metagenomic sequence reads underwent subsequent microbial community profiling using the annotation software MetaPhlAn2 version 2.6.0 with default settings (24). A merged abundance table was created using the MetaPhlAn2 ’utils’ script. The merged abundance table was edited to only include taxa, which MetaPhlAn2 identified up to species level.

The 16S rRNA gene sequences were identified, extracted, and quantified using Metaxa2 version 2.6.0 in paired-end mode with default settings (18). The 16S rRNA gene sequence reads were classified using the mothur version v.1.40.5 “classify. seqs” command with SILVA version 123 as the reference database, with the options cutoff set at 60, probs set as F, and processors set at 8 (25, 26). A custom Unix script was used to create an OTU (operational taxonomic unit) table based on the classifications. Samples with less than 1000 reads mapping to the 16S rRNA gene were filtered out. Metabolic genes were annotated using the HUMAnN2 (HMP Unified Metabolic Analysis Network) pipeline with default settings and with enzyme categories for gene annotation and grouping (27).

We used Bowtie2 for mapping reads to ARG and MGE databases with the options -D 20, -R 3, -N 1, -L 20, and -i S,1,0.50 28). Following annotation, we used SAMtools version 1.4 was to filter and quantify observed ARG and MGE annotations within each sample (29). If both reads mapped to the same gene, we counted the read as 1 match, and if the reads mapped to different genes, both were counted as hits to the respective gene. We used the ResFinder database version 2.1 to search for acquired ARGs and an MGE database, including genes related to or annotated as IS (insertion sequence), ISCR (insertion sequence common region), intI1, int2, istA, istB, qacEdelta, tniA, tniB, tnpA, and Tn916 transposons (6, 30). We searched for plasmid associated genes using the PlasmidFinder database (31).

ARG and MGE count normalization

The Bowtie2 counts for ARGs and MGEs were normalized to the length of each respective gene. Normalized gene counts were then further normalized to the number of bacterial 16S rRNA gene reads obtained from Metaxa2 output, divided by the 16S rRNA gene's length (18). These normalization steps yield an approximation for the number of genes per 16S rRNA sequence for each resistance gene while avoiding bias due to the differential length of the resistance genes. We chose the 16S rRNA gene for normalization instead of library size to account for variation in nonbacterial DNA content in the samples, as explained in Bengtsson-Palme et al. (32). The 16S rRNA gene normalization allows us to approximate the number of resistance genes carried by bacteria, presuming that the copy number of 16S rRNA genes is similar in studied bacterial populations. The 16S rRNA gene-normalized values were used in all downstream analyses. However, we also compared unnormalized data and library size normalized data to validate our choice of normalization, and these yielded similar results to 16S rRNA normalized values (Supplemental Results).

In addition to the 16S rRNA gene, we used the single copy rpoB gene for alternative normalization. HMMER3 and the rpoB Hidden Markow model retrieved from Pfam database were used to count hits (33, 34). The “hmmsearch” command was run with an e^–02 threshold on sequences translated to amino acids in all 6 reading frames. Length normalization was done in the same way as for the 16S rRNA gene. Statistical modeling for rpoB gene counts was done using GLMs and the quasi-Poisson distribution suitable for overdispersed count data (Supplemental Results).

Meta-analysis

We collected metadata for the samples from the National Center for Biotechnology Information Sequence Read Archive or related publications’ supplementary materials (3, 5, 15–17). The metadata used in the meta-analysis are shown in Supplemental Data 3. We also ran all models without the firstborn twin to confirm estimates were not affected by twin pairs, since some of the data sets included twins. We interpreted studies that specified the infants were healthy or that pregnancies were normal as indicating that subjects did not have NEC or preeclampsia. We also analyzed the effects of other clinical parameters, such as maternal diabetes, infant infection status other than NEC, and antibiotic treatment duration, for those subcohorts where data were available. However, they did not significantly improve the models.

We did a preliminary analysis on the long-term effects of being fed formula using data from Bäckhed et al. (5), Baumann-Dudenhoeffer et al. (16), and Shao et al. (17). Supplemental Note 2 lists the European Nucleotide Archive accession numbers for the infants included in the longitudinal analyses from the Bäckhed et al. (5) and Shao et al. (17) studies. The accession numbers of the Baumann-Dudenhoeffer et al. (16) study are in Supplemental Data 3.

GLMs confirmed the independence of the results of the 5 different cohorts and sequencing platform and library preparation method combinations, as including these explanatory variables in the model did not improve the fit (Supplemental Data 4).

We postulated that infants fed formula in addition to fortified human milk would display similar effects as infants fed formula in addition to unfortified human milk. We also assumed that infants fed fortified human milk might differ from infants receiving only human milk. Therefore, we excluded all infants fed fortified human milk that did not receive formula and compared human milk−exclusive diets and formula-containing diets (n = 206). We separately analyzed the effects of human milk supplemented with fortifier compared to an exclusively human milk–based diet. In Bäckhed et al. (5), gestational ages were between 39.3 and 41 weeks. Gestational ages were not reported separately for each infant, and we used an approximate gestational age of 40 weeks for all infants in this cohort. ARG annotation and quantification, as well as Metaxa2 and MetaPhlAn2 community profiling, were conducted for the obtained meta-analysis data set as described above (18, 24).

Statistical analysis: model type selection and model training

All statistical analyses were performed in R version 3.6.1 (R Studio). MetaPhlAn2, Metaxa2, and HUMAnN2 annotation; ARG and MGE mapping results; taxonomy; and metadata files were compiled into individual data objects in phyloseq version 1.28.0 (18, 24, 27, 35). All custom R codes are available from https://github.com/KatariinaParnanen. All figures were plotted with ggplot2 version 3.1.1 (36). The flow of the statistical analysis pipeline is outlined in Figure 2.

An RF analysis was performed using the caret package version 6.0–84 with training control using the “trainControl” command with the option method “cv” and model fitting with “train” command with the option method “rf,” and the importance of the variable was computed using the “varImp” command (37). All variables in the metadata were used in the RF ANOVA importance (for the full list of available variables, see Supplemental Data 1) to model the outcome variable relative abundance of ARGs normalized by gene lengths and 16S rRNA gene counts. Weight, which correlated with both age of the infant and gestational age, and the type of maternal antibiotics, where inadequate numbers of mothers received the same antibiotic, were omitted from the analysis.

We compared different statistical model types to select the best option for modeling the ARG/16S rRNA gene ratio in the individuals from the 46-infant cohort. We examined the data for signs of overdispersion. Because the variance for untransformed data was much larger than the mean indicating overdispersion, we excluded models with normality assumptions. The data had positive continuous values, and the responses seemed to be linear. Therefore, we considered log-transformed normal and gamma distributions as possible distributions for modeling the data. Gamma distributions captured more of the variance than log-transformed normal distributions. Thus, we chose the GLM with gamma distribution and a natural logarithmic link function. The GLM models were fitted with the “glm” command using the option family “gamma” (link = “log”).

Model training for relative ARG abundance was performed with the “step” function, relying on AICs and using chi-squared tests by dropping out 1 variable or interaction term at a time and checking whether the simplified model performed significantly worse than the model with more explanatory factors. Diet was coded as any formula (Formula), any fortifier (Fortifier), or no formula or fortifier (Breast milk). We used the resulting simplified model from model training to analyze the meta-analysis cohorts without the trainset of 46 infants (n = 180) and with the trainset cohort to investigate whether the model estimates change. We also reperformed model selection with all the clinical metadata and all the infants (n = 206) to examine whether formula was retained in the rebuilt model that used all the meta-analysis data.

Cross-validation corroborates that formula increases ARG load in neonates

We next sought to corroborate the impact of diet on the infant gut ARG load by testing our model on subjects from other metagenome cohorts. We included samples from 4 metagenomic studies. The participant flow chart for the subjects is shown in Figure 1. We excluded cohorts matching our literature reach that did not have metadata for gestational age, diet, antibiotic exposure, or infant age, and which included less than 20 infants or had an unbalanced cohort size compared to other cohorts, resulting in the inclusion of 4 cohorts. An additional fifth cohort by Shao et al. (17), a cohort excluded from model testing, was used to investigate the effect on formula feeding in full-term infants during infancy, and the results are described in Supplemental Results (17). The exclusion criteria for each study regarding preexisting maternal and infant health conditions and maternal drug and substance use are available in the original publications. A summary of clinical parameters from the neonates included in the meta-analysis is shown in Table 2 and described in more detail in Supplemental Data 3 and5. Supplemental Data 4 shows the results of variable auto-correlations.

We tested the model built with the train set (ARG load ~16S rRNA counts + gestational age + formula) using the infants from the 4 other included meta-analysis cohorts (Figure 1). Infants with diets having any formula and exclusive human milk were included and infants exposed to fortifier but not formula were excluded (n = 180; Figure 1), following the statistical analysis pipeline outlined in Figure 2. The test confirmed that formula feeding increases the ARG load compared to an exclusively human milk diet. Formula feeding was associated with an approximately 70% increase in the relative abundance of ARGs in the meta-analysis (n = 180; gamma-distributed GLM, fold change for formula-containing diets compared with exclusively human milk–based diet: 1.69; 95% CI: 1.12–2.55; P = 0.013, adjusting for gestational age and 16S rRNA counts; Figure 4A). Formula was associated with an increase of 69% in the ARG load when not adjusting for other factors (n = 180; gamma-distributed GLM, fold change: 1.69; 95% CI: 1.11–2.57; P = 0.015). Gestational age and 16S rRNA counts, which were the other 2 explanatory variables in the model built using the train set, also significantly impacted the ARG load in the test set (n = 180; gamma-distributed GLM, fold change for gestational age/week of gestation: 0.94; 95% CI: 0.91–0.98; P = 0.0013; Figure 4A). The effect of 16S rRNA counts was small but significant (P < 0.01). We also normalized ARGs to the single copy rpoB gene and library size, but the results did not differ from those obtained using 16S rRNA normalized values (Supplemental Results; Supplemental Table 1).

Open in a separate window

FIGURE 4

Effect of formula feeding on neonatal gut resistome in meta-analysis cohorts. (A) Differences in relative ARG sum abundance/16S rRNA gene copy numbers in infants only fed human milk or any formula. (n = 206; n = 180 for other cohorts; n = 26 for infants not fed fortifier in the train set cohort). Formula feeding was associated with an approximately 70% increase in the relative abundance of ARGs in the meta-analysis [n = 180; gamma-distributed GLM, fold change 1.69 (95% CI 1.12, 2.55); P = 0.013, adjusting for gestational age and 16S rRNA counts]. (B) Differences in relative ARG sum abundance/16S rRNA gene copy numbers in infants fed human milk or any fortifier (n = 67). Neither supplementation of the infant's own mother's milk with fortifier or donor milk had a significant effect on the ARG load (gamma-distributed GLM, P = 0.42 for fortifier; P = 0.91 for donor milk, adjusting for 16S rRNA counts and gestational age). (A and B) A regression line is fitted using gamma-distributed GLMs. The y axes have been square root transformed. Dot sizes depict infant ages, with larger dots being older infants. (C) Effect of formula on ARG load in meta-analysis data sets by cohort (n = 206). The boxplot hinges represent 25% and 75% percentiles and centerline the median. Notches are calculated with the formula median ± 1.58 × IQR/sqrt (n). (D) Differences between most abundant ARG and bacterial classes by diet in meta-analysis data sets. (E) Differentially abundant ARGs in infants with different diets using DESeq2 analysis and an adjusted P value cutoff of 0.05 for reported genes. Genes with positive values are more abundant in formula-fed infants, and genes with negative values are more abundant in exclusively human milk–fed infants. (F) Differentially abundant genera in infants with different diets using DESeq2 analysis and an adjusted P value cutoff of 0.05 for reported genera. Genes and genera with positive values are more abundant in formula-fed infants, and genes with negative values are more abundant in exclusively human milk–fed infants. (E and F) The larger the point size, the more abundant the gene is in the samples. Abbreviations: ARG, antibiotic resistance gene; GLM, generalized linear mode; MLSB, macrolide-lincosamide-streptogramin B; rRNA, ribosomal RNA.

Neither supplementation of the infant's own mother's milk with fortifier nor supplementation with donor milk had a significant effect on the ARG load (Figure 4B; n = 67; gamma-distributed GLM, P = 0.42 for fortifier and P = 0.91 for donor milk, adjusting for 16S rRNA counts and gestational age). The meta-analysis cohorts included 1 other study with infants fed fortified mother's milk and donor milk–fed infants (3). Again, these results support the idea that fortifier and donor milk supplementation have less impact than formula feeding on the neonatal gut ARG load. However, the analyses lack the statistical power necessary to draw definitive conclusions (power <0.5 for fortifier and power <0.1 for donor milk).

The median ARG load was higher in the formula-fed infants in all cohorts (Figure 4C), suggesting that formula-fed neonates may exhibit higher intestinal ARG loads, independent of the study design or gestational age. The most prevalent ARG classes differed between the 2 diets, especially for aminoglycoside and macrolide-lincosamide-streptogramin B (MLSB) ARGs, which were more abundant in infants fed formula (gamma-distributed GLM; P < 0.05; Figure 4D). Actinobacteria and Bacteroidia were significantly enriched in human milk–fed infants and Gammaproteobacteria was enriched in formula-fed infants (quasibinomial GLM, P < 0.05, Figure 4D).

The ARGs of formula-fed infants and infants fed exclusively human milk clustered separately (all cohorts, n = 206, including 26 infants not fed fortifier in the trainset; PERMANOVA, Horn-Morisita similarity: R² = 0.011; adjusted P = 0.004). Notably, several ARGs were enriched in formula-fed infants (n = 206; DESeq2: P < 0.05; Figure 4E). The enriched genes included potential ESBLs of the SHV type (P < 0.05; Figure 4E), the mecA gene encoding methicillin resistance in Staphylococcus species, and the MLSB resistance gene ermC encoding erythromycin resistance, typically in Staphylococcus aureus. The SHV, mecA, and ermC genes confer resistance phenotypes, all of which are highly relevant in NICUs (45). Sulphonamide resistance genes of the sul1 type were more abundant in formula-fed infants, indicating the enrichment of class 1 integrons, which contribute to multidrug resistance in bacteria, including Enterobacteriaceae (46).

De novo refinement of the initial model

Because our original cohort was small (n = 46), we reperformed de novo model selection for the gamma-distributed GLMs, modeling ARG loads on all collected variables in both the meta-analysis data set and trainset (n = 206; Supplemental Data 3). Infant and maternal antibiotic use, gestational age, infant age, delivery mode, 16S rRNA gene count, and diet had importance values >0 in our RF analysis, and we included them in the full de novo model. We also added sex, preeclampsia, history of NEC, and the type of antibiotics used, as they might also affect the ARG load. The results of the full de novo model are shown in Supplemental Data 4. The final refined de novo model included formula use, as expected, as well as gestational age, 16S rRNA counts, infant's age, and vancomycin and fortifier use. Formula-fed infants had 57% higher ARG loads (gamma-distributed GLM; n = 206; fold change 1.57; 95% CI: 1.12–2.19; P < 0.009; Supplemental Data 4). The model built using the preterm infant cohort explained approximately 43% of the ARG load variation, and the de novo model explained 74% of the variance (0.77 pseudo-R² calculated using the Cragg-Uhler method).

Consistent with our data set for preterm infants, antibiotic use, which was treated as a binary categorical variable (yes or no), did not improve the model in the meta-analysis data set for neonates (n = 206; chi-squared test: P = 0.96; gamma-distributed GLM, fold change = 0.98; 95% CI: 0.47–2.03; P = 0.96; Supplemental Data 4). However, our study was not designed to investigate the effect of antibiotics specifically, and thus, it should not be concluded that antibiotics do not influence the ARG load. The result merely indicates that for our data set and model, the inclusion of antibiotics did not improve the model fit. Additional results of the possible effect of antibiotics are described in the Supplemental Results. Only 1 of the antibiotics given to the infants, vancomycin, was correlated significantly with the ARG load and was associated with a decreased ARG load in neonates (n = 206; gamma-distributed GLM: fold change 0.53; 95% CI: 0.30–0.92; P = 0.025). Vancomycin is often prescribed to treat infections caused by methicillin-resistant Staphylococcus aureus and Staphylococcus epidermidis, which are often multiresistant.

The authors thank CSC, IT Center for Science Ltd. (Finnish Information Technology Center for Science), for providing the computational resources for the study; Abigail Podany for sample collection; and Lasse Ruokolainen and Otso Ovaskainen for statistical consulting.

The authors’ responsibilities were as follows—KMMP: conducted the formal analysis, validated data, and wrote the original draft; KMMP, JH, RL, JW, SLK, and MPV: conceived the study; KMMP and JH: visualized data; KMMP and MM: determined the methodology; KMMP, RL, JW, and SLK: curated data; KMMP, JW, and CJM: performed the investigation; MM, JH, RS, SR, RL, JW, CJM, SLK, and MPV: reviewed and edited the manuscript; JH, RS, and MPV: supervised the study; SLK: secured resources; MPV: acquired funding and administered the project; and all authors: read and approved the final manuscript.

Author disclosures: KMMP received funding from the MBDP (Microbiology and Biotechnology Doctoral Program) at the University of Helsinki, Finland. SLK received funding from The Gerber Foundation. RS is funded by the Sigrid Juselius Foundation, Finland. All other authors report no conflicts of interest.