Analysis of twelve genomes of the bacterium Kerstersia gyiorum from brown-throated sloths (Bradypus variegatus), the first from a non-human host

Methods

Portions of this text were previously published as part of the first author master’s thesis (https://doi.org/10.11606/D.95.2023.tde-07022024-204134).

Ethic statements

This study did not include any animal experiments, and only non-invasive samples were collected during sampling procedures. All research was performed in accordance with relevant national (see below) and international guidelines and regulations (Mellor, Hunt & Gusset, 2015). Guidelines for handling sloths (Dünner-Oliger & Pastor-Nicolai, 2017) were strictly followed. Animal capture and sample collection (rectal swabs) were conducted between 2015 and 2016 with permission of Instituto Chico Mendes da Conservação da Biodiversidade (ICMBio) and Sistema de Autorização e Informação em Biodiversidade (SISBIO) license number 49627, and approved by the scientific board of Fundação Parque Zoológico de São Paulo. In addition, this research was compliant with Brazilian regulations (Lei Arouca, number 11794, October 8, 2008).

Sample collection and bacterial isolation

Rectal swabs were obtained from wild brown-throated sloths (Bradypus variegatus) for culturable bacterial diversity analysis. Among the cultured bacteria, 12 isolates were identified by 16S rRNA gene sequencing as K. gyiorum in samples of 10 sloths.

Disk diffusion susceptibility test

The Kirby-Bauer disk diffusion test was performed to evaluate the susceptibility of K. gyiorum isolates to antimicrobials (Bauer et al., 1966). Fifteen different antibiotics were tested: amikacin (AMI), amoxicillin plus clavulanic acid (AMC), ampicillin (AMP), cephalexin (CFE), ceftazidime (CAZ), chloramphenicol (CLO), ceftriaxone (CLO), coxycycline (DOX), enrofloxacin (ENO), gentamicin (GEN), ciprofloxacin (CIP), imipenem (IPM), neomycin (NEO), norfloxacin (NOR), and sulfamethoxazole plus trimethoprim (SUT). The inhibition zones of each isolate were measured and compared to the standard defined by Clinical and Laboratory Standards Institute (CLSI, 2018).

Sequencing of 16S ribosomal RNA gene (16S rRNA)

Bacterial growth was promoted by inoculating the samples on Blood Agar with 5% sheep’s blood and Mac Conkey Agar at 36 °C for 48 h. After bacterial growth, the morphology of the colonies was assessed and the different types of colonies were isolated using Tryptone Soy Agar and cryopreserved at −80 °C using Tryptone Soy Broth and 20% Glycerol for further16S rRNA analysis. DNA was extracted from fresh bacterial colonies using Wizard Genomic DNA Purification Kit (Promega, Madison, WI, USA) and used as template for PCR amplification of the V1–V9 region of 16S rRNA gene with primers 27F (5′AGAGTTTGATCMTGGCTCAG 3′) and 1401R (5′CGGTGTGTACAAGACCC 3′) (Nübel et al., 1996). PCR reactions (50 µL) with 30–50 ng of genomic DNA, 2U Taq DNA polymerase (Invitrogen, Waltham, MA, USA) and 0.3 µM of each primer were performed on Veriti 96-Well Fast Thermal Cycler (Applied Biosystems, Waltham, MA, USA) under the following conditions: 3 min at 95 °C, followed by 40 cycles of 2 min at 95 °C, 1 min at 50 °C and 3 min at 72 °C, and a final extension step at 72 °C for 10 min. Aliquots of each reaction were analyzed on agarose gel electrophoresis to confirm the amplicon expected size of ∼1,500 bp. The PCR products were purified using GFX PCR DNA and Gel Band Purification Kit (Merck KgaA, Darmstadt, Germany) and subjected to Sanger sequencing using BigDye Terminator v3.1 Cycle Sequencing Kit (Thermo Fisher Scientific, Waltham, MA, USA) and the above mentioned primers 27F and 1401R or primer 782R (5′ACCAGGGTATCTAATCCTGT 3′) (Chun, 1995) on an ABI PRISM 3130XL genetic analyzer (Thermo Fisher Scientific, Waltham, MA, USA).

The obtained sequences were used to assemble the almost full length 16S RNA gene sequences, which were then searched against the Genbank/NCBI database with BLASTn, allowing taxonomic identification of the isolates as K. gyiorum.

Whole-genome sequencing and assemblies

Purified DNA of the 12 K. gyiorum isolates was sent to the Joint Genome Institute (JGI) for sequencing using Illumina NovaSeq S4 platform producing 2 × 150 pair-end reads. De novo assembly was performed using SPAdes 3.14.1 (Bankevich et al., 2012) and then annotated with the Prokaryotic Genome Annotation Pipeline (PGAP) (Tatusova et al., 2016). Additionally, we retrieved four genomic sequences of K. gyiorum from Genbank public database as of 08 November of 2022. Assembly completeness was assessed with CheckM software (Parks et al., 2015). The basic metrics of all sequences used for this study are shown in Table 1.

Average Nucleotide Identity (ANI) assessment and phylogenomic analysis

The Average Nucleotide identity (ANI) was calculated for all pairs of genome sequences with the aim of determining whether all isolates belong to the same species. We used the python pipeline PyANI v0.2.7 (https://github.com/widdowquinn/pyani) to generate ANIm values using MUMmer for alignment. A heatmap was generated from the pairwise ANI values matrix using the pheatmap package in R v3.5.

For phylogenomic analysis, we used Snippy v 4.6.0 (https://github.com/tseemann/snippy) to generate a core genome alignment. The complete genome sequence of K. gyiorum str. SWMUKG01 was used as a reference (accession number: CP033936). We used Gubbins v2.4.1 (Croucher et al., 2015) to remove recombinant regions prior to the phylogenetic reconstruction. From the core genome alignment, we reconstructed a maximum-likelihood tree using IQ-TREE v1.6.12 with 1,000 bootstrap replicates and based on a GTR+F+I+R8 nucleotide substitution model as predicted by ModelFinder (Nguyen et al., 2015). The ggtree v3.0.4 package was used for visualization and annotation of the phylogenetic tree (Yu et al., 2017).

The population structure of the sloth isolates was further defined using RhierBAPS (Tonkin-Hill et al., 2018) based on a core genome alignment.

Pangenome reconstruction, functional annotation, and gene content variation analysis

The annotated genome assemblies were used as input for Panaroo (Tonkin-Hill et al., 2020). Orthologous genes were grouped using protein sequence identity of at least 95% and query coverage of at least 70%, which determined the core and accessory genes of 16 K. gyiorum genomes. The pangenome was represented on a heatmap using the presence/absence matrix of orthologous genes (OGs). Pan- and core genome curves were generated from the presence/absence matrix of orthologous groups retrieved from Panaroo using the R package vegan.

We used the Clusters of Orthologous Group (COG) database (Galperin et al., 2021) for functional annotation, performing a conserved-domain search with RPS-BLAST against COG domains. The command line for RPS-BLAST was rps-blast–query protein.faa -db COG -out protein_out.out -evalue 0.001 -outfmt 6. COG categories were assigned to each ortholog gene of the core and accessory genome. If a gene was assigned to more than one of the 26 COG categories, it was defined as ‘ambiguous’ category.

We used the presence/absence matrix of accessory genes recovered from Panaroo to perform a principal component analysis (PCA) analysis to assess clustering between isolates from different hosts using ggfortify package in R. BLAST Ring Image Generator (BRIG) v 0.95 (Alikhan et al., 2011) was used to align and compare genomes from different hosts. Variable regions identified by BRIG comparison were extracted and manually checked for accurate annotation using EasyFig version 2.2.273 (Sullivan, Petty & Beatson, 2011).

Detection of mobile genetic elements (MGEs), defense mechanism systems, virulence factors and ARGs

To search for phages in the K. gyiorum genomes we selected two command-line tools: PhiSpy (Akhter, Aziz & Edwards, 2012) and VirSorter2 (Guo et al., 2021). Both tools performed well when evaluated for their ability to predict prophages in bacterial genomes (Roach et al., 2022). We adopted as predictions only those sequences that were predicted to be phage sequences by both tools.

For plasmid prediction, we use PlasmidFinder (Carattoli & Hasman, 2020), whereas Integronfinder2 (Néron et al., 2022) and MGEfinder (Durrant et al., 2020) were used to identify integrons and other MGEs. Defense systems were detected in the genomes of K. gyiorum with DefenseFinder (Tesson et al., 2022) using default settings.

ABRICATE (https://github.com/tseemann/abricate) was used to screen the assembly sequences against the CARD and RESFINDER databases for antimicrobial resistance genes (ARGs) detection and against the Virulence Factor Database (VFDB) for virulence factor identification. We set the percent identity and subject coverage thresholds at 70% and 70%, respectively. The ggtree package was used for plotting a heatmap of presence-absence genes.

General features of K. gyiorum genomes

We isolated and sequenced 12 strains of K. gyiorum obtained from rectal swabs of free-living brown-throated sloths, which were sampled during clinical examination at Fundação Parque Zoológico de São Paulo, between 2014 and 2016. These 12 draft genomes plus four publicly available sequences from the GenBank were analyzed (Table 1). Overall, the genome size of these 16 strains ranged from 3.77 to 3.99 Mbp, with an average length of 3.84 Mbp. The genomic sizes of strains from human origin (3.95–3.99 Mbp) were longer than those strains isolated from sloths (3.77 to 3.81 Mbp) (Fig. S1A). GC content ranged from 62.35 to 62.69%; a small difference is observed between isolates from different hosts (Fig. S1B).

Phylogenetic analysis reveals a structured set of K. gyiorum strains

The ANI values for all pairs of genomes ranged from 98.99 to 99.99%, which means that all 16 genomes belong to K. gyiorum species and are closely related, despite differences in genomic size and GC content (Table 1). A maximum likelihood phylogeny was constructed using the core genome alignment of 16 genomes of K. gyiorum as well as other members of the Alcaligenaceae family: Alcaligenes faecalis strain DMS-30030, Alcaligenes aquatilis strain QD168 and Bordetella bronchiseptica strain NCTC10543, which served as outgroups (Fig. 1A). The phylogenetic analysis showed a host-associated structure in these K. gyiorum genomes, with all isolates from sloths grouping in a single clade, separate from human isolates (Fig. 1B).

The 12 K. gyiorum isolates sequenced in this study were isolated from 10 brown-throated sloth individuals (individuals’ IDs: A, B, C, E, F, G, H, I, J and L). SNP distance and population structure between isolates revealed an important genetic diversity of K. gyiorum in sloths at PEFI. Three clusters (B1, B2 and B3) with at least a thousand SNPs of difference between the clusters were observed (Fig. 1C). Two individuals (B and J) provided two isolates each, and their SNP distance was between zero or two SNPs, which suggests that the intra-host diversity of K. gyiorum may be the result of a single clone, unlike other commensal bacteria, such as Escherichia coli, which maintains multiple lineages within an individual’s gastrointestinal microbiota (Mäklin et al., 2022).

The pangenome of K. gyiorum

We identified 4,199 different orthologous genes in 16 genomes of K. gyiorum, representing its pangenome (Fig. 2A), where 2,974 (70%) genes were shared by all strains (core genome), and 1,225 (30%) genes were present in the accessory fraction (cloud and shell genomes). The pangenome rarefaction curve shows an increasing trend (Fig. 2B). Heaps law allows calculation of whether or not the pangenome is open or closed based on the value of α in the equation n = κN^−α, where n is the number of genes for a given number of genomes (N) and parameters κ and α are calculated in the process of curve fitting (Tettelin et al., 2008). In our analysis, α = 0.91, which is evidence that the K. gyiorum pangenome is open. However, the gene frequency distribution across strains shows a closed pangenome pattern (Fig. 2C), with most genes present in all strains (Domingo-Sananes & McInerney, 2021); this can be explained by the low strain diversity (lowest ANI is 98.99%) that we used in this analysis.

Functional assignation of K. gyiorum pangenome fractions

The COG classifications of each gene of the core and accessory genomes are depicted in Fig. 3 (detailed descriptions are presented in Tables S1 and S2). First, 3,124 (74.4%) of all 4,199 genes from the pangenome were assigned to one COG category, while 451 (10.7%) had an ambiguous COG category annotation (more than one COG category), and 622 (14.8%) lacked a COG annotation (unclassified). The COG categories more frequently present in the pangenome were K (Transcription), E (Amino acid transport and metabolism), P (Inorganic ion transport and metabolism), J (Translation, ribosomal structure and biogenesis), and M (Cell wall/membrane/envelope biogenesis). These five categories accounted for 7.66, 7.31, 5.24, 5.07 and 4.81% of the K. gyiorum pangenome, respectively. The E, K, J, P and I (Lipid transport and metabolism) categories together represent 35.9% of all genes in the core genome (Fig. 3A), whereas categories K, X (Mobilome: prophages, transposons), V (Defense mechanisms) and L (Replication, recombination and repair) are the most frequent in the accessory genome (25.7%). The relative number of genes present in the accessory genome was several-fold higher than in the core genome for the following informative COG assignments: mobilome: prophages, transposons (X, 20.8-fold), Intracellular trafficking, secretion, and vesicular transport (U, 2.44-fold), and Defense mechanisms (V, 2.04-fold).

Specific gene content variation in sloth and human isolates

Principal component analysis (PCA) of the accessory genome shows that gene content variation discriminates K. gyiorum genomes according to host (sloths and humans) (Fig. 4A). We identified 76 genes that were present only in sloth isolates and 54 genes that were present only in human isolates. The functional annotation of these genes revealed different functional assignments (Fig. 4B). While categories I, E, G (Carbohydrate transport and metabolism) and Q (Secondary metabolites biosynthesis, transport and catabolism) are more frequent in genes found only in human isolates, categories X, U (Intracellular trafficking, secretion, and vesicular transport), and V are more frequent in sloth isolates (Tables S3 and S4).

The genes that were found only in the 12 K. gyiorum genomes from sloths were arranged in three genomic regions (RS01–RS03) and in an Inc2 plasmid (Figs. 4C and 5B). RS01 contains poorly characterized genes; RS02 contains genes involved in metabolic function: Energy production and conversion (C) and Carbohydrate transport and metabolism (G), while the RS03 region contains CRISPR-associated genes (cas) that belong to the Defense Mechanisms (V) category (Fig. 4D). The IncN2 plasmid sequence was predicted to be 21,134 bp long (see below).

We found two regions present in all human isolates but absent in sloth isolates: RH01 and RH02, with 18 and 14 genes, respectively (Fig. S2). For both regions, most genes were annotated with COG categories C (Energy production and conversion), E (Amino acid transport and metabolism), G (Carbohydrate transport and metabolism), and I (Lipid transport and metabolism).

Mobile genetic elements (MGEs) and defense systems

IncQ1, IncQ2 and IncN2 replicons were detected with at least 70/70% of identity/coverage sequence thresholds; IncN2 and IncQ2 were present in sloth isolates but absent in human isolates (Fig. 5A). The IncQ1 replicon was detected only in the SWMUKG01 strain isolated from humans, but no plasmids were reported in a previous study (Li et al., 2019). The contig containing the IncN2 replicon, assembled in different ways in the 12 sloth isolates, has direct repeats at the ends, which is evidence that it represents the complete plasmid sequence; it is 21,134 bp long (without the repeats) (Fig. 5B). A BLASTn search of this sequence against the “nt” database returned a match with 70% subject coverage and 95% identity with the pXap41 plasmid found in the bacterium Xanthomonas arboricola pv. pruni (Pothier et al., 2011). This plasmid carries genes encoding type III effectors and helper genes that are absent in the plasmid of K. gyiorum. In contrast, the IncN2 plasmid contained Type IV secretion system genes (virB6, virB5) and other involved in pili formation, plasmid conjugation (traF, G, H, L and J) and segregation (parA and parB) (Fig. 5B).

Twenty-one phage sequences were predicted (Table S5). Nine were predicted by VirSorter2 and 12 by PhiSpy. We considered as phages only sequence regions that were predicted by both tools, resulting in five different phage sequences (Table S6). We gave these the following provisional names: K. gyiorum phage Kgϕ1 through Kgϕ5 (Figs. 5A, 5C). Kgϕ1 (27.7 Kbp) was the only phage detected in sloth isolates (in four isolates). Each of the four human isolates has a different phage sequence (Kgϕ2–Kgϕ5). All phages were predicted to be dsDNA phages. We did not detect any integrases, transposases or insertion sequences (IS) in the K. gyiorum genomes analyzed here.

Microbes can defend themselves against phages and other MGEs using a variety of systems. We detected 14 different defense systems, including cas and restriction-modification (RM) systems (Fig. 5A and Table S6). In all sloth isolates two CRISPR-cas type I operons were predicted, one of them including CRISPR spacers (Fig. 4C). Additionally, other systems were present exclusively in sloth isolates: the defense retro-transcriptase (DRT) system (drt3ab and drt4) and dynamins (LeoA and LeoBC), both involved in anti-viral defense functions. More diverse defense mechanisms were found in human isolates, including RM type I, II and IV; BREX, CBASS, Rst, and Gabija.

Virulence factors and antimicrobial resistance genes

We found 51 putative virulence factors in K. gyiorum genomes. These genes are associated with adhesion, biofilm formation, capsule, flagella, LPS, iron uptake and secretory systems (Fig. 6). Most of these sequences are conserved in isolates from both hosts. The flagellar regulon of K. gyiorum was described by Li et al. (2019) in the SWMUKG01 genome. Here we have identified operons cheARWY, flgCDEFGHIK, flhABC, fliFGIMPORS, and motAB, found in all isolates with the exception of fliR, which was missing in four sloth isolates. Cell surface components like capsule and lipopolysaccharides (LPSs) are essential virulence factors in some Gram-negative bacteria. We identified a highly conserved operon tviBC, whose gene products are involved in the biosynthesis of LPS and capsule in pathogenic bacteria like Salmonella enterica serotype Typhi and Bordetella bronchiseptica (Wear et al., 2022). The wcaJ gene, whose product is related to capsule formation, was only identified in human isolates. On the other hand, five genes, whose products are associated with iron uptake, including the sitABC operon, which is involved in iron uptake in Staphylococcus (Massonet et al., 2006), were found in all K. gyiorum genomes analyzed here. The bcr gene (alcS in Bordetella spp.) was identified only in three human isolates and encodes a permease that maintains appropriate levels of cellular alcaligin levels in Bordetella species (Alcaligenaceae family) (Brickman & Armstrong, 2005).

Phenotypic characterization of antimicrobial resistance in sloth isolates using disk diffusion susceptibility testing revealed susceptibility to almost all 15 antibiotics tested; resistance to norfloxacin (quinolone) was observed in two isolates (Table S7). Genomic inspection using the CARD database detected 20 chromosomal genes putatively encoding multidrug efflux pump systems against fluoroquinolones, tetracycline, phenicols, cephalosporins, as well as other antimicrobials (Fig. 6). These included the ceoAB-opcM encoded efflux pump with resistance to aminoglycoside and fluoroquinolone; mexGHI-opmD encoding resistance to fluoroquinolone and tetracycline, muxABC-opmB efflux pump encoding resistance to aminocoumarin, macrolide, monobactam and tetracycline; and the complex smeABC, encoding resistance to aminoglycoside and cephalosporin.