We investigated 39 pediatric patients with either leukemia (33) or lymphoma (6) from the Atlantic Maritime provinces of Canada, undergoing care at the IWK Health Centre in Nova Scotia, Canada ( Table 1 ). The 33 leukemia patients included 28 acute lymphoblastic leukemia (ALL) patients and five acute myeloid leukemia (AML) patients, and the six lymphoma patients included two hodgkin lymphoma (HL) patients and four non-hodgkin lymphoma (NHL) patients. Specific on subtypes of cancers and treatment protocols for each cancer type are included in Supplemental Material . Patients were asked for stool samples prior to start of chemotherapy and based on convenience while in hospital for treatment, however given the difficulties of collection in this patient population and variable assistance of clinical staff there was no prescribed pattern to sample collection and number of samples and timing varied considerably between patients. Timelines of stool sample collection are in supplemental material ( Supplementary Figure S1 ). A total of 127 stool samples were collected among the 39 patients with a mean of 3 samples per patient (min=1; max=12; and median=2). Pre-chemotherapy samples were collected for 16 patients (8 ALL; 5 AML; 2 HL; 1 NHL) with 12 of these (6 ALL; 5 AML; 1 HL) also having a sample collected during treatment. Among ALL patients 3.4 samples were collected on average from each patient (median=2; min=1; max=12), AML patients had on average 4.6 samples per patient (median=2; min=1; max=11), HL and NHL patients had an average and median of 1.5 samples per patient (min=1; max=2). Stool samples were stored at -20 C for transport and then frozen at -80 C until analysis. Information on antibiotics given, including length of use and number of courses, were collected from each patient. Stool samples were grouped based on the number and duration of courses received prior to the stool sample. All patients receive trimethoprim-sulfamethoxazole throughout treatment prophylactically and so this antibiotic is not included as part of the use and duration. The two main antibiotics used in this patient population are piperacillin-tazobactam (a β-lactam antibiotic) and vancomycin (a glycopeptide antibiotic), although piperacillin-tazobactam is the main β-lactam antibiotic used additional β-lactam antibiotics as well as other antibiotic classes were also given ( Supplementary Table S1 ). Antibiotics were given either intravenously or orally and method of administration is noted in Supplemental Table S1 . Samples were grouped based on multiple-use of: 1) β-lactams, 2) vancomycin, and 3) “any antibiotic” use (AAb). The defined groups within each were: no use (except for trimethoprim-sulfamethoxazole), 1-2 courses, and 3+ courses, of the above 3 groups of antibiotic(s). Additionally, the number of courses of β-lactams, vancomycin, and AAb given were also examined. Samples were also grouped based on duration of 1) β-lactams, and 2) vancomycin. The three defined groups within these two treatments were: no use (except for trimethoprim-sulfamethoxazole), short duration (prescribed ≤10 days), and long (prescribed >10 days) or mixed duration (mix of short and long courses) referred to as long/mixed (LM). Additionally, the number of days of β-lactam and vancomycin antibiotics given were also examined. Timelines of stool sample collection and antibiotic use are in supplemental material ( Supplementary Figure S1 ). The IWK Health Centre Research Ethics Board approved the study and all patients and/or their guardians provided informed consent or assent.

Quality controlled and decontaminated paired-end shotgun metagenomic sequences from a set of 48 matched oral and stool samples from 24 pediatric oncology patients of varying cancer types (see Supplementary Table S2 ) and 5 healthy relatives were assembled using MEGAHIT (kmer values: 21,29,39,59,79,99,119,141; Li et al., 2015) to create a database of contiguous nucleotide sequences. Protein coding sequences were identified in all contiguous sequences greater than 1000 nucleotides in length using PROKKA (version 1.14.6 Seemann, 2014). All protein-coding sequences were compared to the comprehensive antibiotic resistance database (CARD 2021, McArthur et al., 2013) using resistance gene identifier (RGI 4.2.2) with DIAMOND (Buchfink et al., 2015) for alignment with loose hits included. Nucleotide sequences from all identified matches to ARGs were compiled and used to create a bowtie2 pediatric-oncology-ARG-database using bowtie2-build (Langmead and Salzberg, 2012).

DNA was extracted from 127 stool samples using Norgen stool DNA isolation kits as per manufacturer’s protocol. A minimum of one nanogram of extracted DNA was prepared using the Nextera XT library preparation kit (Illumina) following manufacturer’s protocol. Libraries were pooled and sequenced using paired end NextSeq sequencing (150 bp) (Illumina). The Microbiome Helper metagenomic standard operating procedure was followed (https://github.com/LangilleLab/microbiome_helper/wiki/Metagenomics-standard-operating-procedure-v2). Paired forward and reverse reads were run through the kneaddata pipeline for sequence pre-processing. Low quality sequences (reads smaller than 50 base pairs, and with low quality scores PHRED <Q20) were removed using Trimmomatic (Bolger et al., 2014). Bowtie2 (Langmead and Salzberg, 2012) implemented through kneaddata was used to screen out human and PhiX174 contaminant sequences. High-quality, decontaminated forward and reverse paired sequences were concatenated into a single FASTQ file and METAPHLAn3 (Beghini et al., 2021) was used to profile the taxonomic composition of the sequences in each sample. Only taxa present in ≥10% of samples were included in analyses.

The quality-controlled, decontaminated forward and reverse paired sequences from the 127 leukemia and lymphoma samples were mapped to the pediatric-oncology-ARG-database created using bowtie2 (Langmead and Salzberg, 2012). Counts of sequence reads that mapped to each ARG in the database were obtained for each sample using samtools “sort”, “index” and “idxstat”. Mapped read counts were corrected by the number of sequence reads in each sample. While reads were mapped to all ARG sequences identified, only those ≥60% sequence identity were used in downstream analyses. Antibiotic classes were assigned to each ARGs using the CARD database designation, with two exceptions, 1) genes that occurred in an antibiotic class connected with β-lactam drugs were coded as β-lactam antibiotic class genes (i.e., carbapenem, penam, etc.), 2) genes that occurred in multiple antibiotic classes (i.e., penam, fluoroquinolone, glycopeptide), were coded as “multidrug” antibiotic class genes. Counts within samples assigned to the same gene were summed for downstream analysis. Only genes present in ≥5% of samples were used. Genes in four antibiotic classes were selected for closer analysis: β-lactam antibiotic class, glycopeptide antibiotic class, peptide antibiotic class, and multidrug antibiotic class. These classes were specifically selected as the β-lactam antibiotic class and multidrug antibiotic class potentially contains genes for resistance to β-lactam antibiotics, and the glycopeptide antibiotic class, peptide antibiotic class (a parent class to glycopeptide antibiotics), and multidrug antibiotic class potentially contains genes for resistance to vancomycin. All analyses were carried out on gene sequence data, no allele or SNP information was used.

The diversity of species in a community depends on both the number of species observed (species richness) and the relative abundance of species in the community (species evenness). Two communities can have similar numbers of species, but in one community a single species could dominate resulting in differences in evenness between communities. ANOVA and post hoc Tukey’s test was employed to examine changes in richness and evenness in samples with repeated courses (no use, 1-2 courses, 3+ courses of β-lactam, vancomycin, AAb) and increased duration (no use, short, LM of β-lactams and vancomycin) of antibiotics.

Sampling varied from patient to patient and consequently some patients were sampled more than others. Patient samples could differ for phase of chemotherapy, and number and duration of antibiotics. Samples were aggregated over patients and combined into the groups defined above (e.g., number of courses or duration). Note that variation in sampling among groups was not caused by non-random selection or removal of patients. Nevertheless, we wanted to ensure that the groups did not contain a large number of samples originating from just one patient. Hence, we tested for uniform contribution of patients to the samples within a group via the G-test (α<0.05) as implemented in the DescTools R package (Signorell et al., 2022). In each group where the uniform contribution was rejected, we carried out a backward-elimination procedure at the patient level. Specifically, the patient with the largest number of samples was removed (i.e., all samples for that patient were removed only for the group being tested) and the test was repeated for that group. This procedure was repeated until the G-test was no longer significant. In this way patients with a potential overrepresentation to a particular group were identified and removed. While this does avoid the potential for a single patient driving the results, it is expected to lower efficiency in the case where data are missing at random.

Differential abundance analysis of taxa abundance based on number of courses (1-2 vs 3+) of β-lactams, vancomycin and AAb, and duration of courses (short vs LM) of β-lactam and vancomycin were investigated using ANOVA-like differential expression (ALDEx2) (Fernandes et al., 2013). ALDEx2 uses center log ratio transformation and the inter-quantile log ratio (IQLR) was used to determine the features to retain as the denominator of the geometric mean calculation. Using IQLR, accounts for data with systematic variation and centers on the set of taxa that have a variance that is between the lower and upper quartile of variance. The method calculates the expected false discovery rate given biological and sampling variation based on a Wilcoxon Rank Sum test using Benjamini Hochberg (BH) (Benjamini and Hochberg, 1995) (α_BH<0.05).

Differences in mean abundance of the ARGs in samples was examined under three scenarios based on multiple courses of: 1) β-lactam antibiotics, 2) vancomycin antibiotic and 3) “any antibiotic” (AAb) prior to sampling. Mean abundance of each ARG was compared between groups (no use, 1-2 courses, 3+ courses) based on multiple-courses. In addition, abundance of ARGs were compared to number of courses of β-lactam, vancomycin and “any antibiotic” (AAb) prior to sampling. Differences in mean abundance of the ARGs in samples was also examined under two scenarios based on duration of courses of: 1) β-lactam antibiotics, and 2) vancomycin antibiotic prior to sampling. The duration of a course of AAb was not assessed, as the duration was nearly entirely governed by β-lactam antibiotic use. Mean abundance of each ARG was compared between groups (no use, short, LM) based on duration. In addition, abundance of ARGs were compared to number of days of β-lactam, and vancomycin antibiotics prior to sampling. Multivariate linear models using LIMMA (Ritchie et al., 2015) was applied to mean centered, log transformed gene data for all analyses. To account for differences in treatment between cancers and the effect of time in treatment we included cancer diagnosis, and days from start of chemotherapy as confounders. Results were corrected for multiple tests using BH (α_BH<0.05). For the purpose of discovering candidate ARGs and generating novel hypotheses we employ α_BH<0.1 (Mallick et al., 2021).

Genes identified as significant belonging to the β-lactam, glycopeptide, peptide and multidrug antibiotic classes were examined for correlation with taxa present in samples using spearman correlation implemented using corr.test in the psych R package version 2.1.9 (Revelle, 2021). Significance was corrected for multiple tests using BH (Benjamini and Hochberg, 1995) (α_BH<0.05).

Among the 39 pediatric patients, 28 had acute lymphoblastic leukemia (ALL), five had acute myeloid leukemia (AML), two had hodgkin lymphoma (HL) and four had non-hodgkin lymphoma (NHL). A total of 127 stool samples were collected among the 39 patients. The average number of days among all samples between chemotherapy start and stool sample collections was 88 days, (median =50 days; min=-5; max=471). Within cancer types ALL patient samples were on average 101.8 days from start of chemotherapy (median=81; min=-4, max=471), AML patient samples were on average 57 days from start of chemotherapy (median=39; min=-5; max=183), HL patient samples were on average 2 days from start of chemotherapy (median=-1; min=-5; max=12), and NHL patient samples were on average 23 days from start of chemotherapy (median=14; min=-1; max=55). Eight patients (ten stool samples) received no antibiotics prior to collection of stools with remaining stool samples collected following some form of antibiotic. Information on the mean days since start of chemotherapy, days of β-lactam use, and vancomycin antibiotics use, for samples in each group are summarized in Table 1 . Among the 39 patients, 59% experienced one or more gastrointestinal and or blood infection, 26% experienced one or more gastrointestinal infections (predominately C. difficile), 51% experienced one or more blood infection, and 18% both gastrointestinal and blood infections. Among the 51% that experienced blood infections 60% of those patients had at least one infection with a bacterium that was identified within the 127 stool samples (including Streptococcus mitis, Granulicatella adiacens, E. coli, Enterococcus faecalis, Klebsiella oxytoca, Enterobacter cloacae, and Hemophilus influenzae). Due to timing of sample collection however we could not directly examine if there was a link between the presence of these bacteria in the stool and blood infections.

To avoid potential overrepresentation of one or more patients to a particular group, we used the G-test to evaluate the hypothesis (null) that sampling was approximately uniform among patients within each group. Uniform sampling was rejected in several cases, and backward elimination was employed to remove the patient(s) with the largest number of samples until a uniform distribution could not be rejected in those cases. This resulted in the removal of two patients (12 and 11 samples) in the 1-2 courses of β-lactam, one patient (11 samples) in vancomycin no use, one patient (12 samples) in the 1-2 courses of vancomycin, one patient (12 samples) in the 1-2 courses of AAb, and two patients (10 and 9 samples) in the 3+ courses of AAb ( Supplemental Table S3 ). Within the duration groups this resulted in the removal of one patient (12 samples) from the short duration of β-lactam, two patients (10 and 9 samples) from the LM duration of β-lactam, two patient (12 and 9 samples) in short duration of vancomycin, and one patient (10 samples) in the LM duration of vancomycin ( Supplemental Table S3 ). Note that there was no treatment-related selection criteria associated with these patients, nor were they outliers with respect to their treatment. They simply were sampled a disproportionally large number of times throughout chemotherapy. As with any opportunistic sampling design, there will be treatment-independent factors that contribute to greater ease of sampling for some patients (Henry, 1990). Here we analyze the sample trimmed datasets (lists of patient IDs and samples used are provided in Supplemental Table S3 ) as a means to reduce the potential for outlier driven bias even though it entails a cost in statistical efficiency.

Using shotgun metagenome sequence data, 366 taxa were identified in the 127 leukemia and lymphoma samples. This included two Archaea, 348 Bacteria, seven Ascomycota taxa and nine viruses. Mean number of taxa per sample was 42 (min=9, max =115). For species level analyses, only those taxa present in >10% of samples were included, this resulted in 141 taxa (139 Bacteria, 1 Ascomycota, 1 virus) with the mean number of taxa per sample of 38 (min=6, max=89). Among the bacterial taxa, most were members of Firmicutes (83) followed by Bacteroidetes (23), Proteobacteria (16), Actinobacteria (6), and Verrucumicrobia (1) taxa. The distribution of the 141 taxa in the 127 samples are in supplemental material ( Supplementary Figure S2 ).

We assessed changes in richness and evenness in samples with repeated courses (no use, 1-2 courses, 3+ courses of β-lactam, vancomycin, AAb) of antibiotics by using ANOVA and Tukey’s post-hoc test. Richness ( Figure 1A ) differed significantly with repeated courses of β-lactams (ANOVA p=8.2e-5), vancomycin (ANOVA p=0.047) and AAb (ANOVA p=1.3e-4). Figure 1A also shows the individual contrasts with significant differences in richness revealed via post-tests, which are consistent with a general decrease in richness with increased courses of β-lactam and AAb. Although the box plots are consistent with a similar relationship for vancomycin, the maximum pairwise difference (no use vs 3+ courses) was only marginally significance (p=0.052) ( Figure 1A ; Supplementary Table S4 ). In contrast to richness, in no cases did evenness differ significantly with number of courses ( Figure 1B ; Supplementary Table S4 ).

We assessed changes in the differential abundance of phyla and species within the microbial community between 1-2 and 3+ courses of β-lactam, vancomycin, and AAb by using ALDEx2. At the phyla level Bacteroidetes decreased significantly with repeated courses (β-lactam, p= 1.8e-4; vancomycin, p=0.012; AAb,p=4.5e-7; Supplementary Table S5 ) while Proteobacteria increased significantly with repeated courses (β-lactams, p=7.3e-7; vancomycin, p=2.0e-4; AAb, p=5.1e-5; Supplementary Table S5 ). Firmicutes increased significantly with repeated courses of β-lactams (p=0.0053) and AAb (p=0.0035) but not for vancomycin (p=0.1079), whereas Fusobacteria increased significantly with repeated courses of vancomycin (p=0.0486) and AAb (p=0.005) only ( Supplementary Table S5 ).

The species that changed significantly with repeated courses of antibiotics are shown in Table 2 . The relative abundance of 9 taxa differed significantly with repeated courses of β-lactam antibiotics (p_BH<0.05). Four species decreased with repeated courses of β-lactam antibiotics (Eggerthella sp, and 3 species of Firmicutes). Five species increased in relative abundance with repeated courses of β-lactam antibiotics; these included opportunistic Proteobacteria (Klebsiella pneumoniae, and Alcaligenes sp.) and Firmicutes (Enterococcus faecalis, Clostridioides difficile, and Streptococcus). Interestingly, no taxa differed significantly with repeated courses of vancomycin. Examining repeated course of AAb we found that the relative abundance of 18 species differed significantly. Eight species decreased with repeated courses of AAb (Eggerthella sp, Collinsella aerofaciens, Alistipes putredinis, Bacteroides ovatus and 4 species in Firmicutes). Ten species increased with repeated courses of AAb (Rothia mucilaginosa, along with four species in the phylum Firmicutes and five species in the phylum Proteobacteria). Like β-lactam antibiotics, repeated courses of AAb were associated with an increase in opportunistic pathogens.

The overall trend of taxonomic change with repeated course of antibiotics was a decrease in relative abundance in taxa important for gut health and integrity (Clostridium symbiosum, Oscillibacter, Faecalibacterium prausnitzii, Ruminococcus obeum, Eggerthella sp) including those involved in bile acid and short chain fatty acid production. In addition, many of the taxa that increased in abundance with repeated courses (Granulicatella adiacens, Enterococcus faecalis, Streptococcus, Clostridioides difficile and Klebsiella pneumonia) have been implicated in infections in immunocompromised patients.

To examine changes in richness and evenness in samples with increased duration of antibiotics (no use, short, LM of β-lactams and vancomycin) we used ANOVA and Tukey’s test. A significant effect of duration on species richness was only observed for β-lactams, where richness decreased with increasing duration ( Supplementary Table S4 ). There was no relationship between evenness and the duration of either β-lactams or vancomycin treatment.

We investigated the impact of short versus LM duration of antibiotic treatment (β-lactams and vancomycin) on the relative abundance of individual taxa within the microbial community by using ALDEx2. At the phylum level, longer duration of β-lactam was associated with a significant decrease in relative abundance of Bacteroidetes (p=5.2e-4), and increase in Proteobacteria, and Fusobacteria (p=5.7e-4; and p=0.0067; respectively) ( Supplementary Table S5 ). With increased duration of vancomycin only Proteobacteria was significantly changed, showing an increase in relative abundance (p=0.0191; Supplementary Table S5 ). At the species level, we identified no lineages that differed significantly between short and LM duration of β-lactams or vancomycin when we corrected p-values for multiple tests ( Table 3 ).

We identified 264 ARGs spanning 24 antibiotic classes present in at least 1 of the 127 leukemia and lymphoma patient samples ( Supplementary Table S6 ). We found that 254 of these genes occurred in ≥2 samples and 200 of these ARGs occurred in ≥5% of the 127 samples (>6 samples) ( Supplementary Table S6 ). Analyses were performed using the 200 genes (23 antibiotic classes) identified in ≥5% of samples ( Supplementary Tables S6 and S7 ). Most ARGs (69) were connected to resistance in multiple drugs (coded as multidrug), followed by ARGs belonging to the following antibiotic resistance classes: glycopeptide (19), β-lactam (18), tetracycline (18), aminoglycoside (15) and peptide (13) antibiotics. The remaining 49 ARGs were spread across 17 antibiotic classes ( Supplementary Table S7 ). The high occurrence of ARGs in glycopeptide and peptide antibiotic classes, which contain the vancomycin resistance genes, and the β-lactam antibiotic class, which contains resistance genes to β-lactam antibiotics was not surprising, as vancomycin (a glycopeptide antibiotic) and piperacillin-tazobactam (a β-lactam antibiotic) were the most used drugs in our patients (46 and 82% of patients respectively). Based on this, we further investigated the effect of the repeated courses, or duration, of β-lactam and vancomycin antibiotics as well as combined antibiotic use on the 119 resistance genes from the following four antibiotic classes: β-lactam (18 genes), glycopeptide (19 genes), peptide (13 genes), and multidrug (69 genes). These classes were specifically selected as the β-lactam antibiotic class and multidrug antibiotic class potentially contains genes for resistance to β-lactam antibiotics, and the glycopeptide antibiotic class, peptide antibiotic class and multidrug antibiotic class potentially contains genes for resistance to vancomycin.

Because both genes and taxa were impacted by antibiotic use, and it is not possible to directly infer which taxa were contributing ARGs (due to the lack of genomic context, and the potential for LGT between genomes), we used Spearman correlation analysis to investigate associations between the abundance of the 141 species and the 78 ARGs we identified as associated with antibiotic use ( Tables 4 , 5 ). We found significant (p_BH<0.05) associations in 512 of the 10,998 tested comparisons between taxa and ARGs ( Supplementary Table S12 ). We found that over half of these (274) were associated with just 18 taxa across four phyla ( Figure 2A ). The four taxa with the most correlations (Klebsiella pneumonia, Escherichia sp., E. coli, and Enterobacter cloacae) were all from the phylum Proteobacteria ( Figure 2A ). All five species of Proteobacteria we identified were predominately associated with multidrug and peptide genes, and all were opportunistic pathogens. Species in the phyla Bacteroidetes, and Firmicutes were predominately associated with β-lactam and multidrug ARGs in addition to a few glycopeptide genes ( Figure 2A ). The exception to this was the Fimicutes, Clostridioides difficile, Enterococcus faecalis, and Coprobacillus sp., which had associations with peptide genes. Clostridioides difficile in particular had counts of associations with genes in antibiotic classes which were much more similar to lineages in Proteobacteria than to Firmicutes ( Figure 2A ).

Analysis of the significant gene-taxon correlation coefficients showed two distinct clusters of gene-taxon pairs ( Figures 2B, C ). The taxa within one cluster ( Figure 2B ) contained lineages from Actinobacteria, Bacteroidetes and Firmicutes, with most of these taxa having been identified in our previous differential abundance analysis as decreased with repeated courses and duration of antibiotics ( Tables 2 , 3 ). In addition, the ARGs within this cluster were those β-lactam, glycopeptide, and multidrug genes that we previously identified as decreased with repeated courses and duration of antibiotics ( Tables 4 , 5 ). These results suggest that there may exist a unique subset of taxa that decrease jointly with certain ARGs (likely because they encode those ARGs in their genomes) either in response to the sensitivity of those taxa to repeated courses and duration of antibiotics or to the broader ecological changes within the GIT microbiome caused by those antibiotics.

The second cluster ( Figure 2C ) contained lineages from Proteobacteria and C. difficile, Enterococcus faecalis, and Coprobacillus sp., three Firmicutes taxa, two of which increased with repeated courses of antibiotics ( Table 2 ). In addition, taxa in this cluster included associations with peptide genes ( Figure 2A ), and a different set of β-lactam, and multidrug genes than in cluster one ( Figure 2C ). Most of the ARGs in this cluster were genes that we previously identified as increased with repeated courses and duration of antibiotics ( Tables 4 , 5 ) with the only exceptions being vanTG and lsaB. These results suggest that there may exist a different subset of taxa that increase jointly with certain ARGs (likely because they encode them) in response to the ecological changes within the GIT microbiome caused by repeated courses and duration of antibiotics.