Genomic analyses of Burkholderia cenocepacia reveal multiple species with differential host-adaptation to plants and humans

Abstract

Background: Burkholderia cenocepacia is a human opportunistic pathogen causing devastating symptoms in patients suffering from immunodeficiency and cystic fibrosis. Out of the 303 B. cenocepacia strains with available genomes, the large majority were isolated from a clinical context. However, several isolates originate from other environmental sources ranging from aerosols to plant endosphere. Plants can represent reservoirs for human infections as some pathogens can survive and sometimes proliferate in the rhizosphere. We therefore investigated if B. cenocepacia had the same potential.

Results: We selected genome sequences from 31 different strains, representative of the diversity of ecological niches of B. cenocepacia, and conducted comparative genomic analyses in the aim of finding specific niche or host-related genetic determinants. Phylogenetic analyses and whole genome average nucleotide identity suggest that strains, registered as B. cenocepacia, belong to at least two different species. Core-genome analyses show that the clade enriched in environmental isolates lacks multiple key virulence factors, which are conserved in the sister clade where most clinical isolates fall, including the highly virulent ET12 lineage. Similarly, several plant associated genes display an opposite distribution between the two clades. Finally, we suggest that B. cenocepacia underwent a host jump from plants/environment to animals, as supported by the phylogenetic analysis. We eventually propose a name for the new species that lacks several genetic traits involved in human virulence.

Conclusion: Regardless of the method used, our studies resulted in a disunited perspective of the B. cenocepacia species. Strains currently affiliated to this taxon belong to at least two distinct species, one having lost several determining animal virulence factors.

Keywords: Burkholderia cenocepacia; Comparative genomics; Host adaptation; Opportunistic pathogen; PGPR.

Fig. 1

Phylogeny and distribution of host-adaptation genes for 31 B. cenocepacia strains. The evolutionary distances were computed using the Maximum Composite Likelihood method. A total of 1057 conserved core-genes, totaling 1,039,265 positions were used in the final dataset. Branch label colors are indicative of the isolation source of the respective strains. These can either be clinic (red), rhizospheric (green) or environmental (grey). The colored shapes indicate the presence of genetic elements in the genomes of the corresponding strains. Squares correspond to genes that were found to be preferably enriched in clinical (vir.) or environmental (env.) species. From left to right: cable pilus (cblA), 22 kDa adhesion (adhA), Burkholderia cenocepacia epidemic strain marker (BCESM), transcriptional regulator kdgR, bile acid 7-alpha dehydratase (baiE), taurine dehydrogenase (tauX), sulfoacetaldehyde acetyltransferase (xsc), tellurite resistance cluster (telA), low oxygen activated locus (lxa), respiratory nitrate reductase cluster (narIJHGK), nitrate sensor and regulation cluster (narLX), lectin like bacteriocin 88 (llpA), nitrile hydratase cluster (nthAB), phenylacetaldoxime dehydratase (oxd), feruloyl-esterase (faeB), pyrrolnitrin biosynthesis cluster (prn), galacturonate metabolism genes (uxaAB). Circles indicate the presence of the pC3 megaplasmid and the afc cluster. This figure was generated using iTOL [48]

Fig. 2

Whole-genome comparisons of 31 B. cenocepacia strains. The calculations were performed using the Python module PYANI [49]. Two major identity clusters are formed. The bottom cluster consists of B. cenocepacia strains and the second cluster consists of Burkholderia sp. nov. strains. One minor identity cluster is formed by the three outlier strains (Bp9038, CEIB_S5–2, Bp8974) and the last three strains are neither genetically related to B. cenocepacia nor to each other. A double entry heatmap was used to depict the ANI results with ANIm as left entry and ANIb as right entry (a). the dDDH results are depicted on a single heatmap (b). The species demarcation threshold is at ≥95% identity on ≥70% aligned genomic sequence for ANI and at ≥70% identity for dDDH. The exact values and sequence cover ratios are available in Additional file 5: Table S2

Fig. 3

Variations in genomic organization between B. cenocepacia and Burkholderia sp. nov.. The data of 304 genomes presented in Additional file 6: Table S3 was used to represent the differences in genomic organization between B. cenocepacia and Burkholderia sp. nov. strains. Significant levels in variations were determined using Student’s t-test (p < 2.10^− 4, p < 2.10^− 5 for *** and **** respectively)

Fig. 4

Summary diagram of differential adaptation of B. cenocepacia and Burkholderia sp. nov. to different environments. Strains of B. cenocepacia and Burkholderia sp. nov. have been isolated from soils (brown), where they compete with other microbes, plants (green) and animals (red). Burkholderia sp. nov. was repeatedly isolated from soil environment (but also plants, water and aerosols) and possesses several genes improving its fitness in those contexts (green factors). While B. cenocepacia can also thrive in soils, it is often found as an opportunistic pathogen of humans and bears several genes improving its virulence (red factors). Burkholderia sp. nov. can use different plant derivatives (galacturonic acid, xylans, pectin) as carbon sources. It is also proposedly able to synthetize the plant hormone auxin (IAA) through a pathway involving Oxd to convert indole-3-acetaldoxime (IAG) to indole-3-acetonitrile (I3A) which is processed to IAA though the action of Nth. It is also able to produce antibiotics with activity against bacteria (Llpa88) and fungi (pyrrolnitrin). B. cenocepacia strains possess a 22 kDa adhesin which improves its binding to target cells and their invasion. Proposedly, they can also metabolize bile acids, derivatives of cholesterol. In anoxic conditions, B. cenocepacia can survive using its low oxygen activated locus (lxa) and the respiratory nitrate reduction pathway (narIJHG). It also possesses the resistance genes against tellurite, for which the exact functions remain elusive. Source: authors’ design