Transmission of Klebsiella strains and plasmids within and between grey‐headed flying fox colonies

Summary The grey‐headed flying fox (Pteropus poliocephalus) is an endemic Australian fruit bat, known to carry zoonotic pathogens. We recently showed they harbour bacterial pathogen Klebsiella pneumoniae and closely related species in the K. pneumoniae species complex (KpSC); however, the dynamics of KpSC transmission and gene flow within flying fox colonies are poorly understood. High‐resolution genome comparisons of 39 KpSC isolates from grey‐headed flying foxes identified five putative strain transmission clusters (four intra‐ and one inter‐colony). The instance of inter‐colony strain transmission of K. africana was found between two flying fox populations within flying distance, indicating either direct or indirect transmission through a common food/water source. All 11 plasmids identified within the KpSC isolates showed 73% coverage (mean) and ≥95% identity to human‐associated KpSC plasmids, indicating gene flow between human clinical and grey‐headed flying fox isolates. Along with strain transmission, inter‐species horizontal plasmid transmission between K. pneumoniae and Klebsiella africana was also identified within a flying fox colony. Finally, genome‐scale metabolic models were generated to predict and compare substrate usage to previously published KpSC models, from human and environmental sources. These models indicated no distinction on the basis of metabolic capabilities. Instead, metabolic capabilities were consistent with population structure and ST/lineage.


Introduction
The Klebsiella pneumoniae species complex (KpSC) (Wyres et al., 2020) is a group of closely related Klebsiella species which includes Klebsiella pneumoniae, Klebsiella quasipneumoniae subsp. quasipneumoniae, Klebsiella variicola subsp. variicola, Klebsiella quasipneumoniae subsp. similipneumoniae, Klebsiella variicola subsp. tropica, Klebsiella quasivariicola and Klebsiella africana. KpSC are problematic as they function as cosmopolitan opportunistic pathogens (Knittel et al., 1977;Navon-Venezia et al., 2017), responsible for a worrying proportion of community and hospital-acquired infections (Pendleton et al., 2013). Prevalence of multidrug resistance and acquired virulence factors associated with invasive infections is increasing over time (Lam et al., 2021), hence surveillance of potential reservoirs of problematic KpSC sequence types (STs), along with detection of mobile antimicrobial resistance and virulence determinants, is essential for elucidating the potential routes and/or modes of transmission into the human population.
Our previous surveillance study sampled KpSC from four wild grey-headed flying fox colonies in 2017-2018 and represented a snapshot of these populations at the time (McDougall et al., 2021b). Briefly, composite faecal matter was swabbed from drop sheets placed beneath the colonies. Swabs were spatially separated to maximize the likelihood that each represented faecal matter originating from distinct flying fox individuals. KpSC was present in 4/101 samplings from Blackalls Park (flying fox population of 500-2249 on 12/12/2017 and 2500-9999 on 10/4/2018), 18/50 from Camellia Gardens (population of 500-2249), 1/52 from Adelaide (population of 10 000) and 0/52 from Centennial Park wild colonies (population of 10 000-15 999) ( Fig. 1) (Australian Government Department of Agriculture). Additionally, 20 individual captive animals at the Mylor rehabilitation colony were swabbed: KpSC was present in 15/20 animals (total population of 40). In 2021, grey-headed flying fox populations of Adelaide Botanic Park colony were 30 000-36 000 (Boardman, W, pers. comm.). Population estimates for the remaining colonies do not yet have data available. Notably, this was the first non-human sighting of K. africana, a species which has been reported in just two previous studies globally (Rodrigues et al., 2019;Lam et al., 2021). Our initial genomic analyses of these KpSC strains from flying foxes indicated low ST diversity with few STs previously reported among human clinical isolates, as well as low rates of antimicrobial resistance (1.1%) and virulence factors, suggesting that these greyheaded flying fox colonies are not a high-risk environmental reservoir for problematic KpSC isolates. However, this analysis did not explore detailed strain relationships, or the broader genetic content of these KpSC.
Here, we use long-read DNA sequencing to complete the genomes of 13 KpSC from flying foxes and leverage these data alongside the previously reported draft genomes to perform the first high-resolution analysis of KpSC strain and plasmid transmission dynamics within/ between grey-headed flying fox colonies. To better understand the potential for strain and plasmid transfer into the human population, we compare completed plasmid sequences from flying fox isolates to those identified among human-associated KpSC, and utilize the latest genome-scale metabolic modelling approaches to compare the metabolic capabilities of flying fox-derived isolates to KpSC isolates from other sources.

Single nucleotide variant analyses
Single nucleotide variants (SNVs) were identified using RedDog version 1b.10.4 (Edwards et al., 2016). Isolates were grouped by ST, and the corresponding Illumina reads mapped to the respective completed chromosome and plasmid reference sequences (Table S1). As a control, the short-read data from the ST reference isolate were also mapped to its completed assembly. Pairwise SNV distances were calculated for each replicon using the snp2dist script at Figshare (https://doi.org/10.6084/ m9.figshare.16609054) (Vezina et al., 2021). Pairwise chromosomal SNV distances ≤25 were considered indicative of recent strain transmission events, based on previous analyses of KpSC transmission in clinical settings (Gorrie et al., 2017;David et al., 2019;Sherry et al., 2021). An SNV maximum likelihood phylogenetic tree was also constructed from the SNV alignment using RAxML with a GTR substitution model (Vezina et al., 2021).

Plasmid analyses
Completed plasmid sequences were analysed using the mob_typer command from Mob-suite version 3.0.0 (Robertson and Nash, 2018). The closest matching (best query coverage and identity) plasmid sequences from human-derived Klebsiella isolates were identified using NCBI BLASTn search against the public NCBI nucleotide database as of 01/04/21 (Johnson et al., 2008). A MASH distance tree was generated using mash version 2.1.1 (Ondov et al., 2016), along with R packages ape version 5.4-1 (Paradis et al., 2004) and phangorn version 2.5.5 (Schliep, 2011). The tree was inferred using the FastME algorithm (Lefort et al., 2015;Vezina et al., 2021).
Visualizations of pairwise plasmid comparisons were generated via EasyFig version 2.2.5 (Sullivan et al., 2011) and the BLASTn command, with the following BLAST parameters: Min. length -1000. Max. e-value -0.0001, to look at gene content and synteny. Novel metabolic network reconstructions were generated using cobrapy version 0.20.0 (Ebrahim et al., 2013) in a conda environment version 4.9.2 (Anaconda-Software-Distribution, 2016), running Python 3.6.12 (Rossum and Python-Software-Foundation, 1995). The KpSC pan metabolic model (Hawkey et al., 2022) was used as a reference as previously described (Thiele and Palsson, 2010;Norsigian et al., 2020). The specific code used to generate models was adapted from Norsigian et al. (2020) and can be found at Figshare (https://doi. org/10.6084/m9.figshare.16609054) (Vezina et al., 2021). The same approach was performed to update the 22 K. pneumoniae reconstructions described in Norsigian et al. (2019). Model annotations were improved using models iYS1720 (Seif et al., 2018), iML1515 (Monk et al., 2017) and iWFL_1372 (Monk et al., 2013) with a custom script (Vezina et al., 2021). MEMOTE reports (Lieven et al., 2020) were then generated for each model using the report snapshot option. In silico growth capabilities were simulated in minimal media with each of 268 sole carbon, 154 nitrogen, 25 sulfur and 59 phosphorus sources using the flux balance analysis model.optimize method from cobrapy as described previously (Hawkey et al., 2022). All genome accessions and model identifiers are shown in Table S1.
To contextualize the metabolic modelling results with regard to population structure, a core genome phylogeny was inferred from a core gene alignment generated with Panaroo version 1.1.2 (Tonkin-Hill et al., 2020) (options: --mode relaxed -a core --aligner mafft --core_threshold 0.98 -f 0.90 --merge_paralogs). The resulting alignment of 574 339 variant sites represented 3319 genes shared by ≥71/72 isolates. RAxML-NG version 1.1.0 (Kozlov et al., 2019) was used to infer a maximum likelihood phylogeny with 1000 bootstrap replicates (options: --all --model GTR), best of five runs. RAxML was also used to generate an unscaled, daylight phylogenetic tree in the context of plasmid presence (Vezina et al., 2021).

High-resolution SNV analysis identified multiple putative strain transmissions
Among the 13 distinct KpSC STs previously reported from our flying fox isolate collection, five were associated with >1 isolate each (Fig. 2), four of which comprised pairs or groups of isolates that differed by ≤25 SNVs, indicative of recent strain transmissions (Gorrie et al., 2017;David et al., 2019;Sherry et al., 2021). Mapping data can be found in Table S2, while pairwise SNVs can be found in Table S3, ranging from 0 to 715.
There were three occurrences of putative intra-colony transmission at the Camellia gardens colony, including K. pneumoniae ST1017 (n = 2 isolates, three pairwise SNVs), K. africana ST4938 (n = 3 isolates, zero pairwise SNVs) and K. pneumoniae ST5033 (n = 4 isolates, ≤2 pairwise SNVs), and one occurrence at the captive Mylor colony of K. pneumoniae ST5035 (n = 15 of 16 isolates recovered from this colony, ≤3 pairwise SNVs) (Fig. 2). A single putative K. africana ST4938 inter-colony transmission event was observed between the Camellia Gardens and Blackalls Park colonies (n = 2 isolates, 21 pairwise SNVs). Notably, this inter-colony cluster was separated from the Camellia Gardens ST4938 intra-colony cluster and two unclustered isolates by just 36-58 pairwise SNVs, which may be indicative of longer-term circulation within the colony and/or acquisition from a common longterm reservoir (also see ST4938 phylogeny below).
Two K. pneumoniae ST4919 isolates, one from the Adelaide colony and the other from the Camellia Gardens colony (approximately 1150 km apart), were separated by 715 SNVs.

Plasmid diversity and transmission between species
A total of 11 plasmids were identified across the 13 completed genomes. Two isolates (K. africana ST4938, strain FF1003 and K. pneumoniae ST5033, strain FF979) harboured two plasmids each, seven isolates harboured a single plasmid each, while four isolates contained none (Table S2). The plasmids ranged from approximately 5 to 175 kb in length ( Fig. 3 and Table S4). At least seven of the 11 plasmids were predicted to be either selfconjugative or mobilizable and only one contained any antimicrobial resistance genes (K. pneumoniae ST1017, strain FF1009_1, carrying qnrS1 and a variant of dfrA14) and none contained any of the well-characterized KpSC virulence genes (Wyres et al., 2020). Nonetheless, all plasmids demonstrated high identity (≥95%) with generally high coverage (37%-98%; mean 73%) to at least one plasmid from human-derived KpSC isolates (Fig. 3).
A reference-based mapping approach was used to determine the presence or absence of each of our completed plasmids among host isolates of the same ST, which revealed a high level of conservation (≥78% mapping coverage in all but three instances, comprised of three out of seven K. africana ST4938 isolates missing the pFF1003_2 plasmid) (Table S2). At most, plasmids differed among themselves by eight pairwise SNVs (Table S3). Notably, these included the pFF1009_1 AMR plasmid identified in both K. pneumoniae ST1017 isolates, and the 5.5 kb pFF1043_1 plasmid identified among the K. pneumoniae ST4919 isolates collected from two independent colonies (Adelaide and Camellia Gardens). These contained zero and one pairwise SNVs to the completed pFF1009_1 and pFF1043_1 respectively.
Read mapping analyses indicated that a pFF1003_1-like conjugative plasmid was present in nine of the total 39 isolates (23%) (including those from both the Camellia Gardens and Blackalls Park colonies, maximum eight pairwise SNVs). These included representatives of K. africana ST4938, as well as K. pneumoniae ST5037 and ST1412 (Figs 3 and 4). Comparison of the completed sequences for pFF1003_1 (mapping reference from FF1003, a K. africana ST4938 isolate) and pFF1010_1 (from FF1010, a K. pneumoniae ST1412) showed 99% BLASTn coverage and 99% identity, indicative of recent plasmid transmission, though SNV distances were not calculated. However, a comparison between the pFF1003_1 reference and the completed sequence for pFF1023_1 (from FF1023, a K. pneumoniae ST5037) showed only 77% BLASTn coverage and 99% identity. Approximately 9 kb of pFF1023_1 was deleted compared to pFF1003_1 (at least six independent deletion events), including part of the conjugation locus. As such, this plasmid was predicted to be non-conjugative. The pFF1023_1 nan locus [consisting of nanA, nanT, nanE2, nanK, yhcH, a putative Nacetylneuraminic acid porin (GenBank accession: HAT1683918.1), nanC and nanM, flanked by IS3 family transposases] was also inverted compared to the pFF1003_1 reference plasmid (Fig. 4C). pFF1003_1 (97 kb) showed 64% BLASTn coverage and 99% identity to the publicly available 141 kb humanderived KpSC plasmid (pINF116_1). This plasmid was identified from K. pneumoniae ST36 strain INF116, isolated from an Australian human urine sample in 2013 (GenBank accession: CP031793.1). Compared with pINF116, pFF1003_1 lacked many of the amino acid transport genes such as hisQ, hisJ and glnQ, but contained the additional tyrosine recombinase xerC, tRNA(fMet)-specific endonuclease vapC and the putatively intact sialic acid degradation nan locus (Fig. 4C).
KpSC from flying foxes were not distinguished by metabolic capabilities KpSC have specific nutrient and metabolic requirements for growth, which likely influence their host range. To predict whether or not flying fox isolates could potentially colonize humans or other sources and vice versa, we constructed genome-scale metabolic models for the greyheaded flying fox isolates along with KpSC from other sources to identify metabolic commonalities and differences. We then simulated growth phenotypes in silico and overlaid these on the core genome phylogeny of The corresponding heatmap shows the presence of antimicrobial resistance genes in pink, along with identity and coverage to the closest KpSC plasmid isolated from a human source, identified by GenBank's BLASTn. Accession numbers can be found in Table S4. 72 KpSC genomes. All KpSC isolates were predicted to utilize a core set of substrates accounting for 49.6% of carbon, 40% nitrogen, 72.9% of phosphorus and 40% of sulfur substrates tested (Tables S5 and S6). Percentage of total nutrient sources tested which displayed variable usage (non-core) across the KpSC isolates was much lower, consisting of 10.7% carbon, 3.2% nitrogen, 1.7% phosphorus and 0% sulfur sources.
We generated a core gene phylogeny and overlaid the variable substrate usage profiles to examine the relationship between lineage and metabolic substrate usage ( Fig. 5 and Fig. S1). Isolates from flying foxes were distributed throughout the phylogeny and predicted growth phenotypes appeared to cluster by species, then more granularly by ST (Fig. 5), rather than isolate source. Hierarchical clustering of isolates based on predicted growth phenotypes was also consistent with this observation (Fig. S2), with isolates of the same species generally clustering together. One isolate, K. pneumoniae DHQP1002001 appeared as an outgroup to the other K. pneumoniae isolates in the phylogeny. This appeared to be a genetic species hybrid of mostly K. pneumoniae with roughly 1 Mb from a K. quasipneumoniae subsp. quasipneumoniae strain (Fig. S3). Notably, unlike the rest of the K. pneumoniae isolates, this species hybrid was predicted to be unable to utilize ethanolamine and D-glucarate as sole sources of carbon. While K. quasipneumoniae subsp. quasipneumoniae isolates were universally predicted to be unable to utilize the ethanolamine, they were all predicted to utilize tricarballylate and D-glucarate. It is unclear whether this pattern of substrate usage loss exists for Klebsiella hybrid species in general or was limited due to the data used.

Discussion
Our high-resolution genomic comparison of KpSC from grey-headed flying fox colonies across Australia indicates strong evidence for strain transmission both within and between colonies (Fig. 2) because we identified near- Ka refers to K. africana, Kp refers to K. pneumoniae, Kqq refers to K. quasipneumoniae subsp. quasipneumoniae, Kqs refers to K. quasipneumoniae subsp. similipneumoniae, Kqv refers to K. quasivariicola, Kvt refers to K. variicola subsp. tropica, Kvv refers to K. variicola subsp. variicola. Broad ecological group/isolation source is shown by tip colour. Flying fox isolates are shown as light brown strips. Bootstrap values are shown in Fig. S1. Variable predicted growth substrate utilization is shown in the heatmap, coloured by substrate type as indicated. No variable sulfur usage was predicted.
identical strains (very low number of chromosomal SNVs) from spatially separated samples that were likely derived from distinct animals (although it was not possible to associate any given sample with a given individual). While individual flying fox turnover within colonies was not accounted for, the small number of total transmission events was likely due to the rarity of KpSC isolates from flying foxes in general, despite extensive sampling (McDougall et al., 2021b). Five putative strain transmission clusters were identified, four of which involved the wild Camellia Gardens colony in the state of New South Wales. This included three clusters contained within this colony and one transmission cluster between this colony and the wild colony at Blackalls Park. The latter was associated with K. africana ST4938 for which isolates FF1017 (Camellia Gardens) and FF930 (Blackalls Park), were separated by just 21 chromosomal SNVs (Fig. 2) and shared a pFF1003_1-like (98 kb) plasmid with 0 pairwise SNVs, indicating they descended from a very recent common ancestor (although a second plasmid present in FF1017 was absent from FF930). Notably, these colonies are situated 123 km apart (Fig. 1) and were both sampled within a 23-day period. As greyheaded flying foxes have been shown to travel large distances of 50 km (Tidemann and Nelson, 2004) and up to >1000 km in some cases (Roberts et al., 2012), it is possible for individual flying foxes to move between the Blackalls Park and Camellia Gardens colonies contributing bacterial strains to samples derived from both locations and/or causing direct inter-colony transmission. Alternatively, transmission may have occurred indirectly via a common source (food/water etc.). Additionally, our analyses indicated that the two ST4938 clusters, as well as two additional unclustered ST4938 isolates from Camellia Gardens, differed by no more than 58 pairwise SNVs ( Fig. 2 and Table S3). These values are much lower than would be expected for randomly selected KpSC of the same ST (Gorrie et al., 2018;David et al., 2019) and similar to those observed for isolates representing longer-term transmission within hospital networks (Gorrie et al., 2018, David et al., 2019. Hence, these findings are consistent with longer-term transmission in the flying fox population or acquisition from the same long-term reservoir rather than replicate sampling from a single animal, although additional longitudinal sampling of these colonies is required for confirmation. The fifth putative strain transmission cluster was identified at the captive Mylor colony where animals were individually sampled. The colony was located within a rehabilitation centre, which at the time of sampling was a closed population of 30-40 adult grey-headed flying foxes from the local wild Adelaide Botanic Park colony that had been affected by a heat stress event. These unhealthy individuals may have suffered increased susceptibility to opportunistic colonization and dissemination of KpSC isolates, e.g. from other wildlife or domestic animals that were housed in close proximity, which we speculate facilitated the broad-scale dissemination of K. pneumoniae ST5035 in the population [isolated from 15 of 20 sampled individuals (McDougall et al., 2021b)].
There is currently very limited data on the dynamics of bacterial transmission among flying foxes and other bat species; however, our results align with a previous study of viral dynamics (Epstein et al., 2020), which found transmission of Nipah virus in Indian flying foxes within Bangladesh occurred within and between colonies. This study also noted the cycles of Nipah virus prevalence could be influenced by colony demographics, such as individual bat turnover which reduced herd immunity (Epstein et al., 2020). Grey-headed flying foxes in colonies across eastern Australia have an average daily individual turnover of 17.5% as they travel between roosts via a node-based network (Welbergen et al., 2020), indicating high influx and efflux of individuals. This behavioural activity may make grey-headed flying fox colonies susceptible to between-colony transmission of microbes such as KpSC; however, further sampling is required to quantify these transmission events and determine if they represent rare, transient or ongoing phenomena.
In addition to strain transmissions, we identified putative transmission of a 98 kb plasmid, passing the species barrier between K. africana ST4938 and K. pneumoniae ST1412 and ST5037 (Figs 3 and 4). The plasmid variants harboured by FF1003 (K. africana ST4938) and FF1010 (K. pneumoniae ST1412) were highly similar (four pairwise SNVs, 100% coverage and 99% nucleotide identity), supporting a recent transmission event (independent parallel evolution of these replicons is possible, but unlikely). In contrast, the plasmid variant harboured by FF1023 (K. pneumoniae ST5037) showed at least six deletions as well as a structural rearrangement, indicating a plasmid transmission likely occurred at some point in the more distant evolutionary history of these strains (either directly or indirectly via intermediary hosts). In that regard, it is possible that the pFF1003-like plasmid may have been circulating within the Camellia Gardens colony for some time and/or has been introduced from a common reservoir multiple times. Aside from the presence of several selfish genetic elements including a higAB toxinantitoxin system, intact conjugative tra locus (except in the case of pFF1023_1) (Fig. 4C), this pFF1103-like plasmid contained a putatively intact nan sialic acid locus that was conserved in all variants. The nan locus encodes degradation of N-acetylneuraminic acid, a sialic acid (Olson et al., 2013). Utilization and acquisition of Nacetylneuraminic acid improve colonization and provide competitive growth advantages to a variety of pathogens (Chang et al., 2004;Manco et al., 2006;Almagro-Moreno and Boyd, 2009). It is possible that the presence of this locus provided some sort of selective advantage to promote the transmission and longer-term maintenance of the pFF1103-like plasmid in the Camellia Gardens population, although further investigations would be required to confirm this hypothesis.
Plasmid pFF1003 also showed high similarity to a plasmid identified from INF116, a K. pneumoniae clinical isolate causing a urinary tract infection in a patient in a Melbourne hospital in 2013 (Wyres et al., 2019). In fact, our search of publicly available nucleotide sequence data indicated that all completed plasmids from flying foxderived KpSC shared high identity (≥95%) with generally high coverage (37%-98%; mean 73%) (Table S4) to at least one other plasmid identified from human clinical isolates. Hence, our data provide evidence for the flow of genetic material either directly or indirectly between KpSC isolated from flying foxes and human specimens at some point in their evolutionary history. In order for direct transfer to occur these isolate populations would need to coexist, at least transiently, within the same spatial-temporal and ecological niche. Alternatively, indirect gene flow through intermediate plasmid hosts would also explain the presence of these highly similar plasmids. However, we cannot speculate on the directionality, frequency or time frames over which such transfers may occur without larger, contemporaneous sample collections from both flying fox and human populations.
While much more work is required to fully understand the distribution and flow of KpSC strains and genetic material between niches (including detailed genomic analysis of large contemporaneous isolate collections), the results of the preliminary metabolic analyses presented here are consistent with niche overlap (i.e. KpSC from flying foxes did not appear to be metabolically distinct). We hypothesised that there may be metabolic differences between KpSC from flying foxes and those from other sources because occupation of different niches/hosts is likely associated with different nutrient availability, which are impacted by elements such as host diet. We note that our analysis was constrained to a small sample size and that there are likely many metabolic traits which are uncaptured within the metabolic models used in this study, as well as limitations of the current input database. Over time, this approach should become more robust as databases become more complete (https://www.genome.jp/kegg/docs/upd_map.html) (Kanehisa and Goto, 2000). There may also be other genetic differences between flying fox and human strains of KpSC such as K/O loci, previously discussed within the context of flying fox isolates by McDougall et al. (2021b); however, due to a paucity of grey-headed flying fox isolates, we lack the statistical power to discern this.