Metagenome-wide characterization of shared antimicrobial resistance genes in sympatric people and lemurs in rural Madagascar

Sample collection and demographic survey

As a component of a One Health research platform in Ifanadiana District, Madagascar, a household survey was conducted from June to August 2017 in eight communities in roadless areas <5 km from Ranomafana National Park to collect information regarding household member demographics, antibiotic usage, household illness, exposure to wildlife, and previous illness with diarrheal disease. Within these communities, we have documented diverse global health challenges including high prevalence of enteric infections and resistance genes regardless of antibiotic class, and zoonotic human-wildlife linkages (Bublitz et al., 2014; Bublitz et al., 2015; Bodager et al., 2015; Zohdy et al., 2015; Ragazzo et al., 2018). Human participation in the study and collection of survey data were approved of and reviewed by the Emory Internal Review Board (IRB00093812). Before survey administration, informed oral consent was gathered and documented. Household members were also asked to voluntarily submit a fecal sample regardless of history of diarrheal illness. Fecal sample IDs were linked to their corresponding household survey responses, and deidentified for downstream analysis. Fecal samples were collected from captured wild brown mouse lemurs (Microcebus rufus) along footpaths near the villages. The mouse lemurs were trapped using banana-baited Sherman traps (XLR; Sherman Traps Inc., FL, USA), and set overnight at 16:00 and checked at 05:00. One microliter of fresh fecal samples was collected from individual trapped lemurs by using a sterile tongue depressor and transferring the sample into a cryovial filled with approximately 0.8mL RNAlater. The Emory University Institutional Animal Care and Use Committee provided full approval for this research (#3000417) and the field research procedures were approved by Madagascar’s Ministry of Environment, Ecology and Forests (permit nos: 028/17; 083/17; 136/17; 146/17; 164/17).

DNA extraction and sequencing

DNA was extracted from fecal samples using a standard Zymo Inc. bead-beating kit. Whole metagenome shotgun sequencing was performed on the NextSeq 2000 platform using Illumina DNA library preparations. Sequencing produced separate forward and reverse paired-end fastq files, which served as inputs for bioinformatic processing. All mouse lemur metagenomic samples and select human fecal samples from each of the surveyed eight communities were included for metagenomic sequencing and downstream analyses. Human samples were selected through stratified random sampling per village. To be included for selection, households had to include at least one adult, one school-aged child (5–17 years), and one child under 5 years. Eligible households were grouped by their home village and a random three households were selected within each community, sampling without replacement. From the selected households, a sample was chosen for each age group. If only one household member represented an age category, then their sample was selected. If more than one household member was represented by an age group, then a second random sampling was done to choose the representative sample for that age group.

Bioinformatic processing and quality control

An overview of the bioinformatic workflow is shown in Fig. S1. Data quality was assessed with FastQC (v.0.11.9) before and after adaptor trimming and removal of host reads (Andrews et al., 2012). Read quality trimming was conducted using Kneaddata (v0.10.0) with –trimmomatic, which employs Trimmomatic (v0.39-2) (Bolger, Lohse & Usadel, 2014) and Bowtie2 (Langmead & Salzberg, 2012) to remove adaptor reads and reads mapping to the human genome refence GRCh37 (Church et al., 2011), keeping reads at or above Phred 33. Reads from lemur microbiomes were additionally mapped against a draft genome assembly of Microcebus murinus (GCF_000165445.2) (Larsen et al., 2017), selected for its high level of completeness of assembled chromosomes, representation of male and female chromosomes, and shared ancestry to M. rufus. The M. murinus assembly was indexed using bowtie2-build. The forward and reverse paired-end reads from lemur microbiomes were mapped to the indexed assembly using Bowtie2 (v2.5.0) and saved as a SAM file. Unmatched reads (and therefore non-host reads) were subsequently removed using SAMtools’ sort and fastq functions (Danecek et al., 2021). After filtering, only metagenomic sequences with 1 × 10⁶ reads or more were considered for further analysis or assembly. Metagenomes were assembled using SPAdes (v3.15.4) (Nurk et al., 2017) on the trimmed and decontaminated paired-end fastq files using the –meta parameter. Separate forward and reverse fastq files were used as input with the -1 and -2 flags, respectively. The human DNA-scrubbed analysis sequences are available under the Bioproject PRJNA1008138.

Classification of bacterial species and ARGs

Taxonomic composition of the filtered reads was first calculated using MetaPhlAn4 (Blanco-Míguez et al., 2023), with flags –input_type fastq, –unclassified_estimation, and –bowtie2out, to estimate relative abundance of the both classified and unclassified reads which did not match gene markers in the database. This was followed by a second analysis against the Bowtie2 indices to calculate relative abundance of bacterial-associated reads only using the flags –input_type bowtie2out, –t rel_ab, –ignore_eukaryotes, –ignore_archaea. The vOct22 Bowtie2 database available for MetaPhlAn4 was downloaded using the command: metaphlan –install –bowtie2db and used as a reference for taxonomic markers. To calculate the abundance of ARGs, we enumerated the reads per kilobase per million (RPKM) relative to the amount of detected bacterial reads in the sample. We derived this formula from Munk et al. (2022), but accounting for reads. The formula is as follows: ${\displaystyle \frac{\text{Gene reads}}{\left(\text{Length of gene, kilobases}\right)\times \left(\text{Total bacterial reads}\right)}\times 10^9.}$

Filtered reads were first processed through KMA (Clausen, Aarestrup & Lund, 2018) using the AMR Finder Plus nucleotide sequence database to identify ARGs and virulence genes (Feldgarden et al., 2021). ARGs were specifically subset from virulence genes based on classification from the Bacterial Antimicrobial Reference Genes database for sequences related to antibiotic resistance (PRJNA313047). Unique genes were categorized based on annotated gene symbols and unique alleles were categorized based on sequences matching to present NCBI nucleotide reference sequences. Results were included if the detected allele had a template coverage greater than or equal to 60 percent and a query identity greater than or equal to 90 percent, and if they had at least three reads assigned. The results file was then joined with a .mapstat file generated by KMA to quantify the number of reads assigned to each reference sequence. To contextualize the reads relative to bacterial content, the filtered fastq files were also run through Kraken2 (Wood, Lu & Langmead, 2019) using the flags –paired, –report, –classified-out, and –unclassified-out, and referencing the Kraken2 Standard database (26 September 2022) to obtain the number of reads rooted at the bacterial level and the number of unclassified reads, or reads unable to be identified using the database classifications, in the sample. RPKM was calculated for each ARG allele.

Statistical analysis of species and ARG diversity

Final statistical analyses were conducted in R (R Core Team, 2013). Continuous values and counts of discrete data were assessed for normal distribution. The Wilcoxon rank sum test in the stats package (v4.0.4) was used to compare human and lemur metagenomes differences in median total detected reads, proportion of reads mapping to higher order taxa, Shannon indices for bacterial species and ARG allele diversity, RPKM of ARG reads mapping to specific antibiotic classes, and median number of species per sample. Alpha and beta diversity metrics were calculated using the vegan package (Oksanen et al., 2022). Shannon diversity was calculated based on the presence and absence of detected species or alleles. For this system, two measures of alpha diversity were used to help to capture a better understanding of detected species. The Shannon diversity index allows for comparison of both the richness and evenness of the community structure thus relying on both the abundance and the overall number of species, while Chao1 gives a greater weight to low abundance species present in a sample to help predict the likely number of missing species (Kim et al., 2017). Principle component analyses (PCA) were performed on the relative abundances of species and RPKM values of ARG alleles transformed into centered log-ratios to account for the compositional nature of metagenomic data (Gloor, Macklaim & Fernandes, 2016). Subsequently, Aitchison distance was calculated to assess between-sample differences in species/allele diversity. Chao1 and rarefaction statistics were calculated using the iNext package using the sum of the presence of each species or allele detected within human or lemur microbiomes as input (Hsieh, Ma & Chao, 2022).

Differential gene abundance analysis

Differential gene abundance between humans and lemurs for antibiotic resistance genes was conducted on the read counts, summarizing the allele hits to the level of the gene using ALDEx2 (v.1.35.0) (Gloor, Macklaim & Fernandes, 2016). First, the raw read counts were transformed using the command aldex.clr(), with Monte-Carlo sampling set to 128 and the measured denominator set to “all”. To account for both composition and scale in the read counts, uncertainty was added to the model using a gamma value of 0.5. A sensitivity analysis for significance of unique features at various values of gamma was conducted for reads summarized at the gene level and per corresponding antibiotic class associated with resistance (See Fig. S3) (Nixon et al., 2023). To statistically evaluate the transformed abundances, aldex.effect() and aldex.ttest() were used to perform Welch’s T-test and a Wilcoxon rank sum test, and corrected for false discovery using Benjamini–Hochberg corrected p-values (<0.05). Final results were plotted using the aldex.plot() function.

Sequence comparison of antibiotic resistance genes shared between humans and lemurs

AMR Finder Plus (Feldgarden et al., 2021) was used to identify contigs with ARGs. Contigs with the same ARG detected in at least one human and one lemur microbiome were extracted from their sample assembly using bedtools getfasta. These sequences were used to create a custom nucleotide BLAST database (v2.12.0). Each sequence was then queried against the database to identify the pairwise percent identity of each gene compared to other detected genes of the same type. A second BLAST comparison was conducted by extracting the 1,000 base pair regions before and after ARGs commonly present in human and lemur samples, and the pairwise percent identities were quantified. To gain additional insight into the genomic contexts of highly similar ARG-regions, ARG-bearing contigs were extracted from the assemblies and annotated using Bakta (v1.7.0) (Schwengers et al., 2021). For ARG-regions in which commonly found ARGs between pairs of subjects had greater than 90% similarity, the gene synteny of the annotated contigs was inspected and visualized using Gggenes (Wilkins, 2020).

The diversity and abundance of bacterial species and ARGs differ between humans and lemurs

A total of 73 human-derived samples and 15 lemur-derived samples were selected for shotgun metagenome sequencing. Of these samples, 57 human samples had greater than or equal to 1 × 10⁶ total reads after decontamination of human reads. After mapping to the lemur genome assembly, 11 lemur samples had greater than 1 × 10⁶ reads for analysis. The metagenomes of these 68 samples were further characterized for bacterial species abundances and presence of antibiotic resistance genes.

There was a higher total number of reads in lemur samples compared to human samples after decontamination (Fig. 1A) which could not be explained by a large number of small reads present. Although both human and lemur read libraries had average read lengths within an acceptable range for downstream mapping and assembly, the average sequence length for lemurs was higher and tightly ranged (147 to 139 reads) (p < 0.05) (Fig. 1B). For the majority of human and lemur microbiomes, the most abundant taxa identified belong to kingdom Bacteria (Fig. 1C). However, lemur metagenomes noticeably contained higher relative abundances of reads unable to be classified taxonomically. A comparison of unclassifiable reads to taxonomically classified reads did not show significant differences in GC content, Q30 score, average length, or minimum length (Fig. S2). This suggests that much of the diversity of taxa in lemur microbiomes was not represented even in large database collections used for widescale taxonomic composition classification.

Although different in rank of abundance, human and lemur fecal metagenomes include a high abundance of bacteria from the phylum Bacillota (Firmicutes), and human metagenomes are also dominated by Bacteriodota. Several lemur samples were contrastingly dominated by species from Pseudomonadota (Proteobacteria) (Fig. 2A). Humans and lemurs shared 55 of 452 known bacterial species (12.2%) that had at least 0.01% abundance within a single metagenome and occurred in at least 10% of all samples, suggesting there are many rare species detected in the system. The Shannon index between the two groups was higher in humans (P < 0.05) (Fig. 2B). Distinct clustering by sample source was also observed when examining the compositional differences between sample sources. Lemur microbiomes grouped more with other lemur microbiomes without overlapping human microbiomes, and human samples overlapped regardless of the village of residence (Fig. 2C). By Chao1 estimates, when considering detectable species including those in the lowest abundance (<10% of samples) there were fewer species able to be sampled from lemurs (590, 95% confidence interval (CI): 496-729) compared to humans (1091, 95% CI [1057–1144]) (Fig. 2D).

A different pattern emerged when examining the overall composition and abundance of ARGs detected in human and lemur metagenomes. We identified reads mapping to 107 unique ARGs, which comprises 217 unique alleles. Two lemur samples had no reads matching ARGs in our database. Individual microbiomes varied in abundance of ARGs grouped by the class of antibiotic to which they confer resistance, but all human microbiomes carried genes associated with tetracycline and trimethoprim, and 55/57 human metagenomes carried resistance genes to beta-lactam antibiotics. In contrast, there was no one shared antibiotic class among detected ARGs in lemur microbiomes, but all classes seen in human microbiomes were represented in at least one lemur microbiome (Fig. 3A). Lemur microbiomes had no statistical difference in ARG richness compared to human microbiomes (P < 0.05) (Fig. 3B). Still, lemur microbiomes clustered distinctly in their ARG diversity from human microbiomes, however human microbiome ARG profiles from all resident villages overlapped with one another (Fig. 3C). Concordantly, rarefaction estimates, alongside Chao1 calculations suggested that the estimated maximum number of ARG alleles to be sampled are likely similar for lemurs (201, 95% CI [166–264]) and humans (206, 95% CI [160–302]) (Fig. 3D).

Five antibiotic resistance genes were in significantly greater abundance among human microbiomes compared to the lemur microbiomes (Figs. 4A–4C). These genes consisted of dfrF and five separate tetracycline-resistance genes, tet (32), tet (40), tet (W), and tet (Q). The effect size of the difference remained significant for these genes when considering the difference in variance in the data points between the two groups (Fig. 4B). When grouped by associated antibiotic class of resistance, no associated class groups were differentially abundant between humans and lemurs (Fig. S3).

Humans and lemurs share highly conserved integron-associated ARGs

To capture more specific ARG dynamics between humans and lemurs, we quantified and compared assembled ARGs that were detectable in both human and wild lemur metagenomes. A total of 14 ARGs were detected in common between human metagenomes, with all 57 human metagenomes sharing at least one gene with at least one of three lemur metagenomes.

ARGs of the same type were compared between each metagenome containing that gene. Overall, ARGs from different metagenomes were highly similar, with a median nucleotide sequence identity of 99.51 (98.55–99.79) for human-human, 99.67 (99.24-100) for human-lemur pairs, and 100 (100-100) for lemur-lemur comparisons (Fig. 5A). Similarly high levels of sequence identity were found within each ARG (Fig. 5B). The largest ranges of diversity were among pairwise comparisons of tet(O) and tet(Q) genes, respectively (Fig. 5B). We also compared the genetic context surrounding the shared ARGs to determine the similarity of genomic context regardless of ARG sequence conservation. Sequence identity of the 1,000 base-pair flanking regions were also nearly identical between samples. Specifically, the regions around four ARGs (aadA1, dfrA1, qacEdelta1, and sul1) were from seven human metagenomes and one lemur metagenome showed highly similar pairwise sequence comparison of the human-human and human-lemur source pairs (Fig. 5C). Only one gene, lsa(D), had a pairwise identity score less than 100 percent (Fig. 5D).

Because the contexts of these genes were so well conserved, we investigated the gene synteny to understand if there were shared AMR genes in close genetic proximity. ARGs that were co-occurring within assembled contigs included dfrA1 with aadA1, sul1, and qacEdelta. Of the eight samples containing dfrA1, seven contained a conserved aadA1 gene and one sample had a broadly categorized aadA region (Fig. 6). Five of these samples had a sul1 downstream of aadA1, and three contained intI1 (encoding class 1 integrase), including the single lemur sample. Of the seven human microbiomes containing a dfrA1-aadA pairing, residents were from four different villages, with one pairing from the same household. The single lemur sample was collected closest to a village that none of the human residents carrying this cassette were from.