Study population and design
The NHBCS is an ongoing prospective cohort study of over 2000 pregnant women and their offspring recruited from prenatal clinics in New Hampshire, USA. Pregnant women ages 18-45 are recruited from prenatal clinics beginning at approximately 24-28 weeks of gestation and all women must use private well water as discussed previously [52]. Fecal samples were collected prospectively for infants at approximately 6 weeks and 1 year. Institutional review board approval was obtained at Dartmouth with yearly renewal.
Covariate data collection
Extensive covariate data on lifestyle, medical history, and environmental exposures was collected in the NHBCS from participants through medical records, postpartum questionnaires, and telephone interviews. Delivery mode (vaginal or cesarean), feeding mode (exclusively breast fed, exclusively formula fed, or combination of breast and formula fed), prenatal maternal antibiotic usage (yes or no), gestational age (in days), birth weight (in grams), maternal Group B Streptococcus status (yes or no), infant sex (male or female), and other demographic characteristics were collected from telephone questionnaires at 4, 8, and 12 months following birth and delivery medical records. Age of breast milk samples is also used to determine feeding mode and additional information about the classification of breastfeeding exposure. Additional information about feeding method during infant stool sample collection can be found in Additional File 1: Supplementary Methods. Infant antibiotic exposures were only classified if oral, injected, or intravenous antibiotic exposure was indicated during the initial hospitalization. Intrapartum antibiotic exposures administered during labor and delivery were extracted from maternal medical delivery records and were assessed both as a yes/no variable and explored by class of intrapartum antibiotic received. Subjects were grouped according to intrapartum antibiotic exposures using the following categories: no antibiotics; penicillin-like antibiotics only (amoxicillin, penicillin); cephalosporins only (cefazolin, cephalexin); multi-drug classes (two or more antibiotics characterized as penicillin, cephalosporin, vancomycin, clindamycin, and/or gentamicin); or “other” antibiotics such as aminoglycosides, glycopeptides, or lincomycin.
Metagenomics collection and processing pipeline
Infant stool samples were collected at 6-weeks and 12-months maternal postpartum. Stool samples were provided in diapers and stored by guardians in a home freezer (-20 °C) until they were able to return them to the study site. Stool was thawed at 4 °C and aliquoted (range 350-850 mg) into 3ml RNAlater in cryotubes and homogenized before storing at −80 °C. RNAlater stool samples were thawed and DNA was extracted using the Zymo Fecal DNA extraction kit (Cat# D6010, Zymo Research, Irvine, CA), according to the manufacturer’s instructions. For each sample extraction, 400ul RNAlater stool slurry (50–100 mg of stool) was used to isolate DNA. Extractions were performed in batches of multiple samples and included a composite RNAlater stool positive control and a RNAlater negative control. Lysis of stool slurry was performed using 750ul Lysis Buffer in ZR BashingBead™ Lysis Tubes (0.5 mm beads), mixed and then shaken on a Disruptor Genie for 6 min. Eluted DNA was quantified on a Qubit™ fluorometer using the Qubit™ dsDNA BR Assay. Average coefficient of variation of DNA yields (ng/ul) for composite RNAlater stool positive controls was 28%. No DNA was ever detectable in negative control elutions. Concentrations of DNA samples used for metagenomic gene sequencing ranged from 1 ng/ul to 25 ng/ul.
Metagenomic sequencing libraries were prepared at the Marine Biological Laboratory (MBL) in Woods Hole, MA using established methods. These libraries were constructed using Nugen’s Ovation Ultralow V2 protocol. Using a Covaris S220 focused ultrasonicator, DNA samples were sheared to a mean insert size of 400 base pairs.
Some metagenomic samples were processed as paired-end reads (n = 219; 52.1%) and others as single-end reads (n = 201; 47.9%). Generally, sequencing batches of 12 were run for paired-end reads and batches of 16 were run for single-end reads. Paired-end DNA reads were merged into one FASTQ file and all samples were trimmed with KneadData and Trimmomatic with default settings [53] for quality control. Mean (SD) reads after this quality control for single-end reads was 22.3 million (9.6 million) and was 58.9 million (20.8 million) for paired-end reads.
Taxonomic analysis was conducted using both MetaPhlAn2 [35] and PanPhlAn version 1.2.2.5 (10 May 2018) [39]. MetaPhlAn2 [35], an output of the HUMAnN2 [54] version 0.11.2 pipeline, was used to analyze relative abundance of taxa to the species level. MetaPhlAn2 characterizes microbial clades through the use of clade-specific markers. PanPhlAn was used to conduct a strain-level analysis of E. coli. A precompiled pangenome database of E.coli from 2016 [39] was utilized for this analysis. Using PanPhlAn, each sample was mapped to the pangenome to assess presence and absence of genes corresponding to the dominant strain of E. coli in each sample.
ARG markers were quantified using ShortBRED version 0.9.5 [33]. ShortBRED works in two steps. First, it creates a database of antimicrobial resistance gene markers and then uses this set of makers to identify antibiotic resistance in samples. A precompiled list of markers [55] known to confer bacterial antibiotic resistance from the Comprehensive Antibiotic Resistance Database (CARD) [34] was used and the relative abundance of ARG markers in our samples were classified using the “shortbred_quantify” script with default parameters. Outputs from ShortBRED are normalized for average read length, marker length, and sequencing depth and are represented in RPKM. Annotations for the ARGs were adapted from (Sinha et al., 2018), which also used ARG markers from ShortBRED.
MGEs were identified using HUMANn2 [54] with default settings. Briefly, HUMAnN2 uses a ‘tiered search’ approach to first identify known species using reference markers, map these reads to the pangenome of the species, and finally conduct a translated search on all reads that were not classified by known species giving gene family in reads per kilobase. Output gene family files were regrouped by gene ontology (GO), normalized to relative abundance to account for sequencing depth, merged into one file, and then renamed. This was completed using the “humann2_regroup_table”, “humann2_renorm_table”, “humann2_join_tables”, and “humann2_rename_tables” utility scripts available via bioBakery [56].
Baseline metrics of the resistome
The prevalence, mean, and median values for all ARG markers was assessed overall and stratified by 6-week and 1 year samples. Additionally, all values were assessed in comparison to all samples and to only samples that had the gene present. ARGs were considered incident if they were present in at least one sample. A heatmap was created using “pheatmap” [57] to assess the correlation between ARG markers greater than 7,000 RPKM summed across all samples. Relative abundance of these ARGs was then log10-transformed and plotted against the samples. ARGs were clustered using the Canberra distance and samples were clustered using the Euclidean distance based on previously used methods [33].
Covariate selection
Early-life exposures and variables that we hypothesized to be associated with differential resistome composition were chosen based on previous research demonstrating their impact on the resistome and due to their specificity. For our models we used the following covariates: postnatal age of the infant (sampled at approximately 6 weeks or 1 year), gestational age (in weeks), sex (male or female), delivery mode (vaginal or cesarean section), intrapartum antibiotic exposure (yes or no), feeding mode (exclusively breastfed, exclusively formula fed, or mixed fed at the time of sample collection), and infant antibiotic exposure before leaving the hospital (yes or no). Sensitivity analyses assessed sample age in days as a linear variable with the infant as a random effect, intrapartum antibiotic exposures grouped by class of antibiotic prescribed, never versus ever formula fed, and possible joint interactions between intrapartum antibiotic exposure (yes or no) and delivery mode (vaginal or cesarean delivery). Inter-individual differences were analyzed as it often accounts for the largest amount of variation in microbiome studies [20,29,58]. Geographic location has also been associated with differential resistome composition [36], but is controlled primarily through restricting to NHBCS infants.
Statistical analysis
The impacts of covariates on the resistome were assessed primarily using negative binomial regression through quantifying two outcomes: i) the relative abundance of all ARG markers and ii) the presence of unique ARG markers. Additionally, as microbes are hypothesized to be on the causal pathway between many early-life exposures and the resistome, we tested how taxa are correlated with covariates, overall relative abundance of ARGs, and number of unique of ARGs. The outcomes were analyzed across all our samples controlling for sample age (6-week or 1-year) and stratified by the age of sample collection (i.e., cross-sectionally). Negative binomial regression was selected over regression models using the normal or Poisson distribution to avoid overdispersion [19] and because the coefficient can be exponentiated to estimate relative risk. Thus, the interpretation of the negative binomial regression results for assessing the overall number of unique ARGs would be: in comparison to the unexposed group, the exposed group had X times the risk of harboring unique ARGs. Relative risks can be considered statistically significant if they do not include 1, but, if they do cross 1, may be considered meaningfully significant depending on the width of the confidence interval. To assess relative abundance of ARG markers and alpha diversity of species together, Shannon alpha diversity metrics were extracted for species and used as an additional exposure with overall ARG marker relative abundance and number of ARG markers observed as outcomes in sensitivity analyses. To assess the direct correlation between microbial relative abundance and the two overall resistome outcomes, we assessed the Spearman correlation between each CLR-transformed species and phyla relative abundance against each outcome.
Phyloseq [59] objects were created to measure diversity metrics and to make compositional plots for ARGs, taxa, and functional analyses. Alpha (within) sample diversity metrics were calculated using the Shannon and Simpson diversity metrics. For beta (between) sample diversity, genes in RPKM were transformed into compositional data. Using the “microbiome” package [60], a pseudo count of the minimum value divided by two was used in place of any 0s and then the data was CLR-transformed [61]. PCA plots using Euclidean distances were created to visualize results by different covariates. PERMANOVAs were created using the adonis2 function in “vegan” [62] to evaluate between sample diversity.
To further explore similarities and differences in groups with adjustments for all other covariates, MaAsLin2 [38] was used. MaAsLin2 uses a feature reduction technique involving additive boosting of generalized linear models to choose covariates that are most associated with the outcome of interest. Compositional abundance data for each ARG was associated with covariates in MaAsLin2. Deviation from the default parameters included the use of a CLR normalization approach, no standardization of continuous variables, no transformation, and only associations with Benjamini-Hochberg multiple hypothesis correction (q-value) less than 0.01 were considered statistically significant.
To understand whether E. coli strains varied by any covariates in our analysis, we used a variety of visual and statistical tools. Multidimensional scaling using the jaccard distance on the presence/absence matrix of genes with a prevalence of at least 10% was used to visualize sample similarities and differences by covariates. A heat map including only genes that had between 20-80% prevalence was used to cluster samples. Logistic regression with a Benjamini-Hochberg correction (q < 0.01) was used to assess if any genes were differentially abundant by sample age, feeding mode, a combination of delivery mode and intrapartum antibiotic exposure, or total number of reads. A BLAST search of statistically significant genes using KEGG was conducted to reveal differentially prevalent genes.
The goals of the MGE analysis were twofold as we wanted to assess if MGE relative abundance was associated with (1) any covariate or (2) the species-specific contribution of E. coli. Using a custom R script, we collapsed all taxonomic information in the resulting gene families file outputted by HUMANn2 to calculate the relative abundance of each GO term for each infant sample (i.e., the summed gene family relative abundance for each sample equaled 1). Using a CLR-normalization, MaAsLin2 was used to assess which gene families were statistically associated with metadata including sample age, feeding method (ever or never formula fed), intrapartum antibiotic exposure, delivery mode, gestational age, infant sex, and antibiotic use during the initial hospitalization. Statistically significant results (q-value < 0.25) were queried for MGEs. Specifically, the terms “integrase”, “integron”, transpos”, or plasmid” were used for the query as they previously have been used to identify MGE elements from HUMAnN2 results [63]. Since we found these MGEs to only be associated with sample age, we assessed the species specific contribution to each MGE broken down by sample age using the “humann2_barplot” script from HUMAnN2. Bar plots were sorted by the sample sum coefficient of the GO term and grouped by sample age.
Quality control
Number of reads was not added as a covariate in any models where a relative abundance was calculated as all techniques normalized for sample depth (number of reads). However, our assessment of unique ARG markers and genes present in E. coli included log10-transformed number of reads as these presence/absence analyses were impacted by sample depth. Sequencing type and batch effects were assessed through a PCA of the resistome. No evidence of differential effects to the resistome were identified in either (Additional File 1: Figure S6a and Figure S6b) so sequencing type nor batch effect was not considered as a covariate in regression models. A secondary quality control measure was assessed by rerunning a 6-week, single-end sample through the sequencing and downstream analysis pipeline. Upon re-analysis, ShortBRED and MetaPhlAn2 results were nearly identical (Additional File 1: Figure S6c and Figure S6d).