Assessing Transmission of Antimicrobial-Resistant Escherichia coli in Wild Giraffe Contact Networks

We conducted antibiotic susceptibility testing on 765 E. coli isolates from 195 giraffe inhabiting the Ol Pejeta Conservancy (OPC) wildlife reserve in Kenya. Of these isolates, antibiotic resistance was detected in only 2.7% of isolates (n = 21 isolates) from 5.1% of giraffe (n = 10 individuals) (Table 1). Multidrug-resistant isolates (resistant to three or more antibiotic classes [42]) were found in 30% of AMR-positive giraffe (n = 3), with the most common phenotype consisting of resistance to amoxicillin, oxytetracycline, and co-trimoxazole (Table 1). Three isolates from a single giraffe (G670) exhibited cefotaxime MICs above the epidemiological cutoff value (ECOFF) but not the clinical breakpoint, suggesting resistance to cefotaxime may be emerging in wild-type E. coli but is not yet clinically relevant. Similarly, isolate G684₄ exhibited a co-trimoxazole MIC above the ECOFF but not the clinical breakpoint. None of the 765 E. coli isolates exhibited resistance to florfenicol or gentamicin.

Within an individual giraffe, resistant isolates exhibited the same AMR profiles and were typically the same E. coli genomic subtype (mean pairwise genetic similarity based on BOX-PCR: 97.8% ± 1.6%). This suggests that multiple clones of a single E. coli strain were likely collected and screened from individual giraffe. However, between giraffe, genetic similarity was not an indication of matching AMR profiles (mean genetic similarity: 55.6% ± 30.5%).

Given that resistant isolates within an individual giraffe were the same E. coli genomic subtype, a single representative resistant isolate was selected from each AMR-positive giraffe for whole-genome sequencing and de novo assembly (n = 10 isolates). The average read depth per isolate was 28×, with a mean of 248 scaffolds of ≥500 bp after assembly (see Table S1 in the supplemental material). In silico phylotyping revealed that resistant isolates belonged to all four main phylogroups and accessory group E (Table 2). Phylotypes B1 and D accounted for 40% (n = 4) and 30% (n = 3) of the isolates, respectively, while phylotypes A, B2, and E were each represented by one isolate (Table 2). In silico multilocus sequence type (MLST) analysis identified seven E. coli sequence types (STs), including ST1670, which accounted for all three of the phylotype D isolates (Table 2).

A total of 10 different acquired ARGs were identified, with an average of 4 ARGs per isolate (range: 0 to 8 ARGs [Table 2]). Genes conferring resistance to sulfonamides had the highest prevalence (sul1 and sul2; n = 9 isolates), followed by genes conferring resistance to tetracyclines [tet(A) and tet(B); n = 8], aminoglycosides (strA and strB; n = 7), and beta-lactams (bla_TEM-1; n = 5). In general, ARGs identified in each isolate reflected the observed phenotypic AMR profiles (Tables 1 and 2). However, we did not identify any acquired ARGs in isolate G670₁ that could be responsible for observed resistance to amoxicillin or the elevated cefotaxime MIC. Additionally, isolate G366₂ carried genes that typically confer resistance to oxytetracycline [tet(A)] and co-trimoxazole (sul1, sul2, and dfrA7) but did not exhibit phenotypic resistance to either antibiotic. Finally, although the presence of sul2 in isolate G684₄ could account for a co-trimoxazole MIC above the ECOFF, the same gene was also identified in two other isolates that did not exhibit elevated co-trimoxazole MICs (G688₁ and G363₁) (Tables 1 and 2).

Sequence alignments of ARG-carrying scaffolds identified six scaffold groups where pairwise scaffold similarity was ≥99% for at least 90% of the shorter scaffold’s length (Table 2). Two of these scaffold groups were characterized by the presence of a sul2-strAB gene cluster (scaffold groups 2 and 5), while scaffold group 1 carried the genes catA1, dfrA7, and sul1 (Table 2). The remaining three scaffold groups (3, 4, and 6) carried one ARG each: tet(A), tet(B), and bla_TEM-1, respectively (Table 2). Genes associated with mobile genetic elements (MGEs), such as those encoding transposases, resolvases, integrases, and plasmid replicons, were found in 71% of all ARG-carrying scaffolds (n = 17) across 50% of isolates (n = 5) (Table 2). Of particular note was the presence of the class 1 integron-integrase gene, intI1, on all members of scaffold group 1, in addition to a unique scaffold carrying ARGs tet(A) and dfrA5 (Table 2). Interestingly, the group 1 scaffold from isolate G366₂ also carried the phage-related integrase gene, intS, suggesting the presence of a phage. Six of the 17 ARG-carrying scaffolds were identified as plasmids across five isolates, including the IncQ1 plasmid in three sul2-strAB isolates and the IncK plasmid in one bla_TEM-1 isolate (Table 2). All members of the bla_TEM-1-positive scaffold group also carried the gene tnpR, which encodes a transposon gamma-delta resolvase typically found in Tn3 family transposons.

We investigated the importance of host-level social and ecological factors for predicting antimicrobial-resistant E. coli in giraffe using univariate and multivariate Firth bias-reduced logistic regressions (Tables 3 and 4). In both univariate and multivariate models, social degree (i.e., number of social contacts) was the single most important predictor of AMR risk, accounting for more than 90% of Akaike’s information criterion corrected for small sample size (AICc) weight in univariate models and included in all three best-fit multivariate models (Tables 3 and 4). Specifically, giraffe with lower social degree were significantly more likely to carry at least one resistant E. coli isolate (P value ≤ 0.001). Age class was the second most important predictor in univariate models, with 5% of AICc weight (Table 3). Neonates were significantly more likely to carry at least one resistant E. coli isolate than juveniles, subadults, and adults (Tables 3 and 5). However, in multivariate models, only one of the three best-fit models contained age (Table 4). This result was somewhat surprising given that 71% of all neonates (5/7 giraffe) carried at least one resistant isolate and neonates accounted for 50% of all AMR-positive giraffe (Table 5). One possible explanation for the minimal importance of age in regression models is that other factors, such as social degree, may have reflected additional information about age variation not captured by age classes. For example, among giraffe with known ages (i.e., neonates and juveniles), social degree was highly correlated with age (Pearson’s r = 0.80; Pearson’s r = 0.93 if restricted to <150 days of age) (Fig. 1), and younger individuals were more likely to carry resistant E. coli (Wilcoxon rank sum test: W = 3; P value < 0.0001).

AMR risk was also related to small home range size and low spatial degree in univariate models (Table 3), but neither factor was consistently present or significant in the best-fit multivariate models (Table 4). Although social and spatial betweenness accounted for less than 1% of AICc weight in univariate model comparisons (Table 3), both factors were included in all three best-fit multivariate models (Table 4). Specifically, giraffe AMR risk increased with increasing social betweenness and decreasing spatial betweenness (Table 4). However, the values of betweenness regression coefficients were small (Table 4), suggesting that the effect of social and spatial betweenness on AMR risk was at most minor compared to other factors.

We next assessed the relevance of giraffe contact networks for the transmission of ARGs between giraffe. Plots of both social and spatial networks suggest there was relatively little clustering of AMR-positive giraffe (Fig. 2). Network k-tests and path tests (43) were conducted on social and spatial networks for each scaffold group identified in at least two AMR-positive giraffe. All tests failed to reject the null hypothesis that ARG-positive giraffe were distributed randomly in the networks (Bonferroni adjusted P values > 0.99).

This study was conducted at Ol Pejeta Conservancy (OPC), a wildlife reserve located in Laikipia, Kenya (0˚N, 36˚56′E) that integrates commercial cattle ranching with wildlife conservation. The reserve is bisected by the Ewaso Ng’iro river, with the western side home to wide-ranging OPC cattle herds and the eastern side featuring small clusters of cattle and proximity to local villages and their livestock. All giraffe within OPC at the time of the study (n = 212) were individually recognized based on unique spot patterns on their necks. Age class (neonate, <3 months; juvenile, 3 months to 1.5 years; subadult, 1.5 to 4 years; and adult, >4 years) for each giraffe was established according to physical attributes and age-associated behaviors (21, 78, 79). Approximate birth dates were known for all neonates and juveniles. The population exhibited a 50:50 sex ratio, although the sex of two neonates born at the end of the observation period was not determined. At the end of the study period, the giraffe population of OPC consisted of 160 adults, 20 subadults, 21 juveniles, and 11 neonates.

Behavioral observations were carried out in OPC from 21 January to 2 August 2011 as described elsewhere (21). Briefly, giraffe were located by driving daily routes through different regions of OPC. Giraffe social group membership was determined by proximity to other group members and/or movement of individuals in a common direction (35, 89, 90). In total, there were 1,089 sightings of giraffe groups during the study period, with each giraffe observed an average of 31.1 ± 7.6 (mean ± standard deviation [SD]) times (approximately once per week).

This research was approved by Kenya’s National Council for Science and Technology (permit NCST/RRI/12/1/MAS/147) and the UC Davis Institutional Animal Care and Use Committee (protocol no. 15887).

We constructed two giraffe contact networks based on (i) social associations and (ii) spatial overlap of home ranges, with the same set of individuals included in both networks (n = 193). For the social network, patterns of association were established from group membership, with pairs of giraffe linked if they were observed in the same social group at least once. Network connections were weighted according to association strength (AS), which was calculated as the number of sightings where a giraffe pair was observed in the same group divided by the total number of times they were seen together or apart (21). Because giraffe mothers tend to isolate their calves from other giraffe for the first 1 to 3 weeks postpartum (78, 79), two neonates born at the end of field observations were linked only to their mothers.

The spatial network was constructed based on the extent of home range overlap between individuals (21). Network connections were weighted according to the proportion of overlap. Home range boundaries of each giraffe were mapped using a fixed-kernel utilization distribution of Global Positioning System (GPS) coordinates recorded during sightings. Each giraffe’s core home range was calculated using a 75% kernel density estimation (91). Average home range size ranged from 16.9 ± 13.4 km² for neonates to 95.7 ± 3.3 km² for adult males.

For both social and spatial networks, we calculated two standard measures of network connectivity for each individual: weighted degree and weighted betweenness (92). Degree is defined as the number of individuals to which the focal individual is connected (93). Weighted degree (here social/spatial degree) is the extension of degree for weighted networks and accounts for both the number of linked individuals and the weight of those links (92). Betweenness measures the extent a focal individual falls on the shortest paths between other pairs of individuals in the network (93). For weighted networks, the shortest paths used to calculate betweenness are based on the sum of connection weights (92). Previous work has shown that both metrics positively correlate with microbial diversity and bacterial subtype sharing (21, 22, 94), as well as pathogen infection risk (95, 96), suggesting that individuals with high degree or betweenness may have elevated opportunities for exposure to antibiotic-resistant bacteria compared to those with low network connectivity.

Detailed methods relating to fecal sample collection, E. coli isolation, and genetic analysis are described elsewhere (21). Briefly, fecal samples from 195 giraffe were collected between 10 August and 11 September 2011 (after the completion of behavioral observations). Collected samples were streaked onto E. coli selective CHROMagar EC agar (CHROMagar, Paris, France) and incubated overnight at 37˚C. Up to six E. coli colonies (range, 2 to 6; median, 4) were randomly selected from each incubated sample, subcultured, and frozen for transport to the United States. To determine the genetic profile of each isolate, densitometric curves were generated using the banding patterns from BOX-PCR and gel electrophoresis (97–99). Genetic similarity of each isolate to all others was subsequently calculated through pairwise comparisons of densitometric curves using the Pearson product-moment correlation coefficient. Two isolates were interpreted as the same E. coli subtype if they were at least 90% similar (21).

All E. coli isolates (n = 765) were tested for sensitivity to six antibiotics: amoxicillin (beta-lactam: penicillin), cefotaxime (beta-lactam: third-generation cephalosporin), florfenicol (phenicol), gentamicin (aminoglycoside), oxytetracycline (tetracycline), and co-trimoxazole (sulfonamide/dihydrofolate reductase [DHFR] inhibitor). These antibiotics were selected based on reports of antibiotic use in Kenya and recent research in eastern and southern Africa demonstrating acquired resistance in E. coli from local livestock and wild animals (8, 27, 48, 49, 65, 68). In addition, amoxicillin and oxytetracycline were the primary antimicrobial treatments used by OPC in cattle and wild animals (G. Omondi, unpublished data).

Sensitivity testing was performed using broth microdilution in a two-step process following methods established by the Clinical and Laboratory Standards Institute (CLSI) (100) and European Committee on Antimicrobial Susceptibility Testing (EUCAST) (101). First, isolates were tested against a single concentration of each antibiotic based on the epidemiological cutoff values (ECOFFs) determined for Enterobacteriaceae by EUCAST (102). ECOFFs represent the MICs that separate wild-type bacterial populations from strains with acquired resistance mechanisms (102). While these values do not provide information on the clinical relevance of resistance per se, they are useful for detection of AMR emergence in wild-type populations (102). Thus, an isolate was described as “tentatively resistant” if bacterial growth was observed in the presence of the antibiotic at the ECOFF concentration. In step two, all isolates tentatively classified as resistant from step one were tested for sensitivity across a range of 10 antibiotic concentrations using standard broth microdilution methods described by Wiegand et al. (103). Any isolate with a MIC above the established EUCAST clinical breakpoint was considered resistant (104). E. coli ATCC 25922 was used as a quality control strain in both steps of the sensitivity screen, as recommended by the CLSI and EUCAST (100, 101).

DNA was extracted from one resistant isolate per giraffe using the Qiagen DNeasy Blood & Tissue kit (Qiagen, Hilden, Germany). DNA libraries were prepared and multiplexed with the NexteraXT library preparation kit (Illumina, San Diego, CA) and subjected to paired-end whole-genome sequencing on the Illumina MiSeq platform with 300-bp read lengths. Raw sequences were trimmed and quality filtered using Trimmomatic v 0.33 (105), and sequence quality was assessed with FastQC v 0.11.7 (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Average read depth and genome coverage were determined by mapping reads to the E. coli K-12 strain MG1655 chromosome (NCBI reference sequence NC_000913.3) using the Burrows-Wheeler Aligner v 0.7.17 (106) and SAMtools v 1.5 (107). De novo genome assemblies of high-quality reads were performed with SPAdes v 3.10.0 (108) following the assembly selection procedure of Ahlstrom et al. (47).

Prokka v 0.7.17 (109) was used to identify open reading frames (ORFs) on assembly scaffolds of ≥500 bp. In silico phylotyping of resistant isolates was based on the Clermont scheme (110) and performed using the web-based service ClermonTyper (http://clermontyping.iame-research.center/) (111). In silico MLST was used to identify isolate sequence types and performed using the Center of Genomic Epidemiology web-based MLST 2.0 method (112) (https://cge.cbs.dtu.dk/services/MLST/). Acquired resistance genes were identified by aligning ORFs to the ResFinder database (downloaded 22 April 2018) (113) using BLASTN (114) with an E value threshold of 1e−10, query identity of ≥99%, and query coverage of ≥85%. The similarity between ARG-carrying scaffolds from different isolates was assessed by pairwise BLASTN alignments. Scaffolds were assigned to the same scaffold cluster if their alignment exhibited ≥99% similarity across ≥90% of the length of the shorter scaffold sequence. To determine the genetic context of identified acquired ARGs (e.g., on a plasmid, within a transposon, etc.), we identified genes associated with mobile genetic elements on each ARG-carrying scaffold using the INTEGRALL (115), ISfinder (116), PlasmidFinder (117), and NCBI nucleotide (https://www.ncbi.nlm.nih.gov/nucleotide/) databases.

To investigate whether any host traits influenced the risk of carrying antibiotic-resistant E. coli, we performed univariable and multivariable logistic regression models using Firth's penalized likelihood approach implemented in the R package logistf (118–120). This method corrects for biases in maximum likelihood estimation due to rare events and small sample sizes while still producing finite and consistent parameter estimates (118, 119). The dependent variable (0/1) was whether an individual carried at least one E. coli isolate resistant to at least one antibiotic (i.e., AMR risk). Covariates included age class, sex, home range size, river side (west versus east), social and spatial degree, and social and spatial betweenness. River side was included due to different management strategies for cattle on each side of the river, as well as because previous research has shown that E. coli population structure is in part influenced by the river (121). In a separate univariable logistic regression, we also tested for a relationship between AMR risk and the total number of isolates collected per giraffe to confirm that more isolates did not increase the likelihood of detecting antibiotic-resistant E. coli (P = 0.49). Due to the nonindependence of network data, univariable and multivariable regression P values were calculated via permutation methods where the dependent variable was randomized relative to the covariates (10,000 total permutations) (122, 123). P values were defined as the proportion of random permutations that yielded a coefficient as extreme as the observed value. Covariates with P values of <0.2 in univariable regressions were included in multivariable models. We used forward selection for inclusion of covariates to determine the best-fit multivariate models. To account for potential collinearity among covariates, we calculated generalized variance inflation factor (GVIF) values for each multivariable model and rejected models exhibiting GVIFs of >2.0 (124). GVIFs are recommended over traditional variance inflation factor values when any covariate in the model has more than 1 degree of freedom (125). Candidate models were compared using Akaike’s information criterion corrected for small sample size (AICc) (126). Models with a ΔAICc of <2.0 were considered equally parsimonious (126). All analyses were conducted in R v 3.4.2.

To determine whether giraffe social or spatial association networks represent potential transmission pathways for antimicrobial-resistant E. coli between individuals, we used the network k-test procedure implemented in R (43). This method involves calculating the mean number of infected individuals occurring within one step of an infected individual in a network (i.e., mean infected degree), called the k-statistic. The network location of infected individuals is then randomized (node-label swap) and the k-statistic is recalculated after each data permutation (10,000 total permutations). The P value is calculated by comparing the observed k-statistic to the distribution of null k-statistics. If the mean number of infected nodes within one step is significantly greater than expected if cases were randomly distributed in the network, this suggests that the pattern of infection cases may have resulted from transmission via network links. Due to the large number of connections between individuals in both the social and spatial networks (many of which had low weight and probably represent transient interactions) and because the k-test does not account for link weight, pairs of giraffe were connected only if the link weight (i.e., AS or home range overlap) was greater than the median weight of all network pairs. For weighted networks, we also performed an extension of the k-test procedure that uses the average weighted inverse path-length between each infected node and the nearest other infected node (which we term the “path test”), with the assumption that infected nodes will be closer together (shorter path lengths) than random expectations if the network represents potential transmission pathways. For this study, we defined an infected case as a giraffe harboring an E. coli isolate with both an ARG-positive scaffold and the corresponding resistance phenotype. Separate k-tests and path tests were conducted on both social and spatial networks for each scaffold group that was identified in at least two AMR-positive giraffe. To correct for multiple hypothesis testing on the same association network, we applied Bonferroni adjustments to all k-test and path test P values within each network.