Core/pan resistome and main ARG structuring factors
Samples from the twelve infants in the cohort were selected in order to represent the temporal distribution of microbiome taxonomic profiles during the first year of life (Fig. S1). Shotgun sequencing resulted in an average number of DNA sequence reads per sample of 121,121,684 (±32,477,513 s.d.). A total of 1682 ARGs were identified using the read-based approach with ShortBRED, with a mean carriage of 140±27 ARGs in the 12 infants, and a mean percentage of 0.27 (range 0.06-1.30; Table S1) reads mapped to CARD. Of the 260 uniquely identified ARGs, 47 were found in all 12 infants, and 47 were unique to a single child (Fig. 1a, Table S2). The core resistome was dominated by efflux systems, particularly systems associated with the TolC outer membrane channel. For example, the AcrAB, EmrAB and mdtABC efflux complexes all depend on TolC for function and are typically associated with E. coli. On the other hand, ARGs unique to a single individual functioned predominantly through drug inactivation. Both in terms of the number of ARGs identified and the overall mapping rate, efflux systems were the most prevalent type of resistance mechanism (Fig. 1b and c). However, while the mapping rate was five times higher for efflux systems relative to the second most prevalent category (drug inactivation) the corresponding ratio for the total number of unique ARGs was only 1.2, demonstrating that many ARGs are associated with non-dominant taxa.
The most important structuring factor was the individual (R2=0.21, p<0.001, PERMANOVA) when looking at ARG carriage in the full cohort. However, there was also a highly significant overall temporal progression in ARG composition, both when looking at the data set as a whole (R2=0.05, p<0.001, Fig. 1d), and when looking at each of the twelve infants individually (p<0.01 for all tests). These results were also corroborated with random forest models. A regression model with sampling day since birth as response variable resulted in an R2 of 0.66 (p<0.001), while a classification model with child identity as the response variable resulted in a classification success rate of 80.5% (p<001). Interestingly, the classification model was unsuccessful in discriminating the individuals comprising the twin pair, ID10 and 11, with a combined error rate of 54%. Together, these results indicate uniqueness of ARG carriage amidst overall general temporal trends in gut resistome composition.
ARG carriage declines over the first year of life
The overall main drug class to which resistance was the most prevalent was beta-lactam antibiotics (Fig. 2). This was also the case at the individual level in eleven of the twelve infants (Fig. S2). Tetracycline was the predominant drug resistance class for ID2, while beta-lactam was only the fifth most prevalent drug class in this individual. When splitting beta-lactam into sub-classes (carbapenem, penam, penem, cephalosporin, cephamycin and monobactam) the most abundant resistance class in seven of the infants was fluoroquinolone. Genes for tetracycline resistance were dominant in four infants, while cephalosporine resistance was the most prevalent class in ID12 (Fig. S3). Resistance to sulphonamides and glycopeptides was rare, and not detected at all in several infants.
In addition to large temporal fluctuations, the infants showed highly individual temporal trajectories with respect to carriage of ARGs conferring resistance to various classes of antibiotics (Fig. 3). Relative abundance patterns of ARGs were quite stable (Fig. S4), but with a tendency towards expansion of tetracycline, lincosamide and streptogramin resistance genes at later time points. The overall percentage of reads mapped to the CARD database showed a highly significant linear downward trend over time (R2=0.22, p<0.001, generalized additive model) (Fig. 4a). ARG diversity also showed a significant, albeit slightly non-linear, negative time trend (R2=0.18, p<0.001) (Fig. 4b). Significantly reduced abundances over time were observed in ARGs representing 13 out of the 18 main antibiotic drug classes (Table 1), while resistance to the two rare classes, glycopeptides and sulphonamides, showed the opposite trend, and streptogramin resistance genes did not show any significant time trend. ARGs within three of the six main resistance mechanisms (antibiotic efflux, target alteration and reduced drug permeability) also decreased significantly in abundance over time (Table 2).
ARG carriage correlates with Proteobacteria abundance and E. coli in particular
The infant gut microbiome is highly dynamic over the first year of life, with variability in population abundances on multiple levels of taxonomy. Leveraging our longitudinal sampling effort allowed us to evaluate contemporaneous correlations between CARD mapping rates to the population dynamics of specific taxa. The percentage of sequence reads mapped to CARD generally correlated strongly with relative abundances of Proteobacteria (mean Pearson R2=0.81; Fig. 5), with significant contemporaneous correlation observed in all children except for ID12 (Table 3). In this child, the CARD mapping rate did not correlate strongly with any particular phylum, and the highest correlation was with Actinobacteriota (R2=0.57, p=0.053). ID12 had by far the lowest average CARD mapping rate (Table S1) of the 12 children, as well as the lowest mean relative abundance of Proteobacteria (1.82%, Table S3) (Fig. 5), and excluding this infant resulted in an increased mean correlation of 0.87 between the CARD mapping rate and relative abundance of Proteobacteria.
In the full cohort, there was a significant decrease in the population of Proteobacteria over time (Spearman’s rho=-0.45, p<0.001). E. coli was the dominant species within that phylum (Fig. 5; Table S4) in the majority of the children, and the mean contemporaneous correlation between CARD mapping rate and E. coli relative abundance was R2=0.69. These correlations were statistically significant for nine infants, marginally significant for ID8 and 12, and non-significant for ID6 (Table 3) where hardly any E. coli was observed (Fig. 5). In ID6, the CARD mapping efficiency tracked the dynamics of Proteobacteria, but this could be attributed mainly to the Klebsiella species K. pneumoniae (R2=0.99) and K. quasipneumoniae (R2=0.97). In ID12, unlike the other infants, we were unable to identify specific taxa with which the ARG mapping rate correlated strongly.
Although at differing degrees of prevalence, E. coli was observed in all twelve infants, mostly with some consistency (Fig. 5), and as demonstrated above this species displayed high ARG carriage. By using matching criteria of at least five non-zero observations of gene and taxon, and a cutoff at 0.97 Pearson correlation, we assigned ARGs to E. coli by matching gene dynamics with the dynamics of the E. coli populations in 11 of the 12 infants. Since E. coli was only observed in ID6 at the four final time points, at a mean relative abundance of 0.05%, we were not able to associate any ARGs with E. coli using this method for this infant. Instead, in ID6, the high contemporaneous correlation with the dynamics of Proteobacteria ARG mapping rates was found to be driven to a large degree by Klebsiella spp. (Fig. S5). For the 11 other infants, an average of 26 ARGs were linked to E. coli (range 5-41; Fig. S5), with 43 unique ARGs identified in the set (Table S5).
The majority of ARGs were annotated on Proteobacterial MAGs
In addition to using a read based approach and contemporaneous correlation to link ARGs with taxa, we also assembled reads into contigs and binned contigs into MAGs. Since contamination from closely related genomes could lead to the artefactual occurrence of multiple ARG copies on a MAG, each unique ARG was counted only once per genome. CARD annotation, with default rgi settings, of all MAGs in all the children combined resulted in 7150 hits. Interestingly, the distribution of percent amino acid identity of hits showed a bimodal distribution, with a roughly normally distributed group of low identity hits centered around 35-40% amino acid sequence identity (Fig. S6). When focusing on this group by only looking at hits with less than 50% identity (3736 hits) we found that 86.7% of these hits were annotated as conferring resistance to glycopeptide antibiotics. More specifically, antibiotic target-altering enzymes conferring vancomycin resistance. We suspect that these hits are false positives that are actually representing widely distributed enzymes with various functions in producing cell wall precursors.
In order to remove false positive hits we applied filtering on 80% fragment coverage and amino acid identity, which left 1989 ARGs that were used for subsequent analysis, with a per infant mean of 166±53. On the phylum level, 80.5 % of these ARGs were found on MAGs classified as Proteobacteria, followed by Firmicutes (15.0%) and Actinobacteria and Bacteroidetes (both at 2.0%). On the species level, 17 out of the 20 top ARG carriers were classified as Proteobacteria (Table 4), with E. coli having by far the highest carriage at 28.0% of all annotations. The three species in the top 20, that were not members of the Proteobacteria, were Staphylococcus aureus (4.0%) Enterococcus faecalis (2.3%) and Staphylococcus epidermis (1.1%). In all twelve infants, E. coli was the dominant carrier of ARGs.
The majority of ARGs assigned to E. coli by contemporaneous correlation are found on E. coli MAGs
Since E. coli had the highest carriage of ARG by a large margin, we focused our analysis on E. coli in order to compare ARG assignment by contemporaneous correlation with MAG annotation. First, the estimated proportion of reads mapped to E. coli MAGs was highly correlated with E. coli population dynamics estimated with MetaPhlAn, with a mean Pearson’s R2 of 0.95, demonstrating good agreement between the two approaches. In fact, the mean correlation was pulled down by ID6 (R2=0.68), without which, the mean correlation was 0.97. Next, we compared ARGs assigned to E. coli by contemporaneous correlation in an infant with ARGs annotated on the corresponding E. coli MAGs. E. coli MAGs were generally of good quality (mean completeness, contamination and strain heterogeneity of 90.7, 11.1 and 73.4%, respectively; Table S6). MAG completeness was reduced in ID10 and ID11, the twins, likely due to the presence of a second Escherichia MAG that could not be reliably classified to the species level, and that likely confounded the binning process. In the other 10 children, a single Escherichia MAG was found and classified as species coli (mean completeness 95.8%). The mean number of ARGs annotated to the E. coli MAGs was 44 (range 32-51), with 13 genes found in all 12 infants and 7 in only a single individual (Fig. S7). In addition, the majority of the ARGs assigned to E. coli by contemporaneous correlation were also found on the corresponding MAG in each child (mean 72.6%), with the exception of ID6, where no genes were assigned by correlation.
A total of 29 unique ARGs were assigned to E. coli by contemporaneous correlation in at least one child without being found on the corresponding E. coli MAG in the child in which the correlation was observed (Table S7). Most of these were found in only one or a few infants, but an AmpC beta-lactamase gene (ARO:3004290) conferring resistance to penicillin-like and cephalosporin-class antibiotics was assigned to E. coli by correlation in ID1, 2, 3, 4, 8, 10 and 11. The fact that this ARG was identified in seven infants, and not found on the corresponding E. coli chromosomes, suggests that it may be associated with a plasmid, or some other mobile element that was not readily assembled on a contig or successfully binned. As another example, in six of the infants we associated acrA (ARO:3004042) with E. coli using contemporaneous correlation, while this gene was not annotated on the corresponding E. coli MAGs. This gene is normally associated with an efflux system in Enterobacter cloacae.
Virulence factor carriage declines over the first year of life, and correlates with Proteobacteria, and Escherichia coli in particular
We also investigated the prevalence of virulence factors (VFs) in the metagenomes of the 12 infants. We identified a total of 4487 genes encoding VFs in all 12 infants combined, and 742 unique VF genes in the full dataset, with a mean mapping rate to VFDB of 0.48% (range 0.02-2.65%). 73 of these genes were shared between all 12 infants, while 131 were unique to a single child (Fig. S8). Like for ARGs, the strongest structuring factor for VF carriage was identity (R2=0.22, p<0.001, PERMANOVA), while time of sampling also had a highly significant overall effect (R2=0.07, p<0.001). As with ARGs, the proportion of reads mapping to the VF database showed a significant general negative time trend over the first year of life (p<0.001, GAM; Fig. S9a). This trend was significant in 10 of the 12 VF categories observed, with stress-survival and motility being the only exceptions (Table S8). Unlike for ARGs, we did not see a significant trend in VF diversity (Fig. S9b). The read mapping rate to VFDB showed a significant contemporaneous correlation to the CARD mapping rate in all 12 infants (mean R2=0.87, Table 5), as well as a strong correlation (mean R2=0.82) with the relative abundance of Proteobacteria (Fig. S10, Table 5). As for ARGs the exception was ID12, where the correlation was negative and non-significant, while the relationship was significant in the other infants (Table 5). Overall VF carriage also correlated strongly with E. coli relative abundance (mean R2=0.75; Fig. S10), with a significant relationship for all infants except ID6 and 12 (Table 5). Furthermore, excluding ID6, a mean of 61 individual VF genes (range 10-137) showed strong (R2>0.97) contemporaneous correlation with E. coli population dynamics (Fig. S11).
Contigs with co-localized ARGs and MGEs decrease in abundance over time
The potential of ARGs to be mobilized and transferred from one taxa to another was investigated by annotating the contig assemblies for both ARG and MGE markers in order to identify physical linkages between the two gene categories. In total, 2726 ARGs were annotated on 2002 contigs, and 729 of these contigs also carried one or more putative MGEs. Kaszab et al. [41] used a universal cutoff of 10 kb separation between an ARG and an MGE in order for an ARG to be considered to have a high probability of being mobilizable. This cutoff was considered conservative relative to those used by other investigators [42, 43]. For the 5 MGE categories defined by MobileOG we modelled the cumulative percentage of ARGs linked with MGEs as a function of the distance between ARGs and the nearest MGE (Fig. 6a). The total proportion of annotated ARGs located within 150 kbp of an MGE was 41%, with more than half of these found in relatively close proximity (<10 kbp; Table 6). Considering all MGE-linked ARGs, irrespective of MGE category or the distance between ARGs and MGEs, normalized abundances showed a significant overall negative trend over time (Fig. 6b). The same time trend was seen for contigs with ARGs and MGEs in close proximity (<10 kbp) (Fig. 6c). Significant negative time-trends were also observed when considering plasmids and phages separately (Table 7). In contrast, IGEs showed a significant positive time-trend. However, it should be noted that the categories Plasmid and Phage (557 and 310 contigs, respectively) showed a substantial overlap (190 contigs). This is due to both dual MGE categories for certain MobileOG genes, like ratA and hokG (104 contigs), and annotation of both plasmid and phage genes on the same contig (163 contigs). The MGE category IS co-localized with ~1% of the total ARGs and was excluded from further analysis.
Of the 491 contigs with at least one ARG in close proximity of an MGE, 79.0% were binned to a MAG. On the phylum level, most of these contigs were binned to MAGs classified as Proteobacteria (66.8%), with 7.4% and 3.9% in bins classified as Firmicutes and Bacteroides, respectively. Escherichia coli carried by far the largest proportion of mobile ARGs with 30.5%, followed by Enterobacter hormaechei A and Citrobacter freundii (both ~5%, Table S9).
MGE categories are differentially linked with different categories of antibiotics and resistance mechanism
The MGE category to which ARGs are linked might influence how frequently these genes are transferred, and to which taxa. Furthermore, different resistance mechanisms could preferentially associate with certain MGE categories. Looking at the distribution of ARG resistance mechanisms by MGE category, and only considering ARGs located within 10kb of an MGE, we did not see significant differences in the distribution of resistance mechanisms between the total ARG pool and the total pool of MGE-linked ARGs (Fig. 7). However, the resistance mechanism profiles within the various MGE categories did differ significantly compared to the total pool of MGE-linked ARGs in some cases (Table S10). For example, in contrast to the other MGE categories, ARGs linked with IGEs did not predominantly encode efflux system components, and the dominant resistance mechanisms of IGE-linked ARGs was drug inactivation. The frequency of linkage between ARGs conferring resistance through target protection and both categories of transposon (IGE and CE) was more than fivefold that of the full distribution of MGE-linked ARGs (p<0.05 for both comparisons). Plasmids and phages had a significantly higher proportion of linkages with efflux pumps (p<0.01), and a lower proportion linked to drug inactivation genes (p<0.01), compared with the full distribution. We also observed substantial variation among individuals (Fig. S12), with for example linkage between the resistance mechanism ‘reduced permeability’ and the MGE categories plasmids and phages only observed in one infant.
Since efflux systems frequently confer resistance to multiple classes of antibiotics, these were excluded from analysis of potential linkages between MGEs and ARGs within different drug classes. With respect to the main drug classes, the distributions of all ARGs mapped to contigs and MGE-linked ARGs showed a few notable difference (Fig. S13, Table S11). ARGs conferring resistance to rifampicin were found on a number of contigs, but in no instance were these genes linked with any category of MGE. In addition, the proportion of ARGs conferring resistance to aminoglycoside and tetracycline were significantly lower in MGE-linked ARGs (p<0.006), while resistance to beta-lactams was higher (p<0.004). Noteworthy differences (Table S12) between MGE categories and the full distribution of MGE-linked ARGs included the lack of tetracycline resistance genes on plasmids, and a high proportion of these genes linked with IGEs compared to the total MGE-linked ARGs (p< 0.05). On the other hand, plasmids were more frequent carriers of beta-lactam resistances genes (p<0.05).