Burkholderia species in human infections in Mexico: Identification of B. cepacia, B. contaminans, B. multivorans, B. vietnamiensis,B. pseudomallei and a new Burkholderia species

Background Burkholderia sensu stricto is comprised mainly of opportunistic pathogens. This group is widely distributed in the environment but is especially important in clinical settings. In Mexico, few species have been correctly identified among patients, most often B. cepacia is described. Methodology/Principal findings In this study, approximately 90 strains identified as B. cepacia with the VITEK2 system were isolated from two medical centers in Mexico City and analyzed by MLSA, BOX-PCR and genome analysis. The initial identification of B. cepacia was confirmed for many strains, but B. contaminans, B. multivorans and B. vietnamiensis were also identified among clinical strains for the first time in hospitals in Mexico. Additionally, the presence of B. pseudomallei was confirmed, and a novel species within the B. cepacia complex was documented. Several strains misidentified as B. cepacia actually belong to the genera Pseudomonas, Stenotrophomonas and Providencia. Conclusions/Significance The presence of different Burkholderia species in Mexico was confirmed. Correct identification of Burkholderia species is important to provide accurate treatment for immunosuppressed patients.


Introduction
second group comprised 26 strains isolated in Hospital Infantil de México Federico Gómez (HIMFG) from children with pneumonia and initially identified as B. cepacia [12]. These strains were included in this study to identify them at the species level ( Table 1). The type species B. cepacia ATCC 25416 T , B. contaminans LMG 23361 T , B. multivorans LMG 13010 T , B. vietnamiensis TVV75 T and B. cenocepacia J2315 T and reference strain B. cenocepacia TAtl-371 were included in the analyses. The strains were regularly kept in LB Petri dishes at room temperature and preserved with 35% glycerol at -70˚C.

16S rRNA and atpD gene amplification and phylogenetic analysis
Total DNA was isolated by growing the strains in 5 mL LB and incubating overnight at 30˚C with shaking (120 rpm). DNA was obtained using the cetyl trimethyl ammonium bromide (CTAB) method according to Moore and Dowhan [23]. The 16S rRNA gene was amplified with primers 27F/1492R [24]. The atpD gene was amplified according to Baldwin et al. [25]. The PCR fragments were sequenced at Macrogen Inc. (https://www.macrogen.com) using the same primers for amplification. The sequences were analyzed and assembled with Chroma-sPro 2.

Polymerase chain reaction amplification of the BOX dispersed-repeated motif (BOX-PCR) analysis
All strains were analyzed with BOX-PCR using the BOXA1R primer [28]. Cycling conditions were 95˚C for 5 min and then 35 cycles of 95˚C for 1 min, 63˚C for 1 min and 72˚C for 6 min and a final elongation cycle of 10 min at 72˚C. The BOX patterns were observed in 1% agarose gels.

Whole-genome sequencing and analysis
Selected strains were grown in 40 mL LB and incubated overnight at 30˚C. The total DNA was isolated using the method of Moore and Dowhan [23]. The genome sequence was obtained at Novogene (https://en.novogene.com/) using Illumina Platform PE150 with a 350-bp insert DNA library. For the genome assembly, the quality control analysis was carried out with FastQC to 0.11.9 (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/), checking for the presence of adapters and length of reads. The clipping and filtering of the information was performed with Trimmomatiq 0.39 [29], keeping only those sequences with phred quality greater than 28. All samples were assembled with an optimized method of on de Bruijn graphs and automatic adjustment of the K value with the SPAdes 3.14 program [30]. Some of the assemblies underwent a cleaning process to remove spurious sequences using Blast 2.9.0 [31] comparing to a Burkholderia genome bank. Metrics such as N50 and misassembles were obtained using QUAST 5.0.2 [32]. Annotation was performed using the standard operating procedure at DOE-JGI Microbial Genome Annotation Pipeline (MGAP v.4). Average Nucleotide Identity (ANI) calculations were obtained with JSpeciesWS [33], using ANIm (based on the MUMmer algorithm). The digital DNA-DNA hybridization (dDDH) was calculated in silico with the Genome-to-Genome Distance Calculator (GGDC 2.1) using the BLAST method [34]. Results were based on recommended formula 2 (identities/HSP length), which is independent of genome length and appropriate for incomplete genomes.

Phylogenomic analysis
A phylogenomic analysis was carried out to determine the relatedness of Burkholderia strains with Burkholderia species using the program up-to-date bacterial core gene (UBCG), which is based in the phylogenetic analysis of 92 core genes, that were previously selected from the analysis of 1429 species available at EzBioCloud (https://www.ezbiocloud.net/), representing 28 phyla, including Burkholderia species [35]. Single-nucleotide polymorphism (SNPs) was analyzed with strain 294H and compared to other B. pseudomallei strains. The analysis was performed with maximum parsimony based on core SNPs with Parsnp, a component of the Harvest 1.3 software (https://github.com/marbl/harvest).

Multilocus sequence typing (MLST) analysis
MLST allele sequences were analyzed by BLAST searches using the sequenced genomes with the Bcc and B. pseudomallei PubMLST Databases [36].

Antibiotic susceptibility
All strains were tested with the VITEK2 system for antibiotic susceptibility. , including those from HGMGG and HIMFG, and the type species of Paraburkholderia unamae used as an external species, showed that the strains were not clearly grouped with any of the species from the Bcc, despite the high similarity found with the 16S rRNA gene analyses, except for a group from HIMFG that was close to B. pseudomallei/B. mallei strains (S1 Fig). Since the phylogenetic analysis with 16S rRNA genes has low taxonomic resolution, which is in line with several reports [37,38,39], a phylogenetic analysis was performed with atpD gene, using the same number of strains. The atpD gene was used because it is one of the genes listed in the PubMLST for the taxonomical analysis of Bcc (https://pubmlst.org/ organisms/burkholderia-cepacia-complex). Phylogenetic analysis of the atpD gene gave a better resolution, showing that the strains from HGMGG formed a large cluster, and the strains from HIMFG clustered in few groups but none clearly close to any Burkholderia species (S2 Fig). The concatenation of both 16S rRNA and atpD genes to seek a better resolution and identification of the strains, resulted in a large group containing all strains from HGMGG (Fig 1). Strains 184D and 662D from HIMFG were close to B. vietnamiensis, with similarity of 99.28 and 99.93%, respectively, and 99.2% similar to each other. Other HIMFG strains, such as 500H and 501H, grouped with B. cenocepacia but the 16S rRNA showed more similarity to B. contaminans (99.85%). The strains identified as B. pseudomallei/B. mallei were intermingled with these species (Fig 1). The positioning of the rest of the strains from HIMFG was not clear, showing no association to any particular Burkholderia species.

BOX-PCR analysis
Since the resolution in the analysis of 16S rRNA and atpD genes was poor in general, the BOX profiles were analyzed visually to determine clonality and to define possible Burkholderia subgroups among the isolates. Twelve BOX patterns were found among the HGMGG and HIMG strains (Table 1); as expected, none of these patterns was shared among the strains from both hospitals. Strains collected from HGMGG displayed four different patterns with pattern I present in the largest number of strains (31 strains). Strains from HIMFG displayed eight different patterns, with pattern VIII present in 10 strains. At least one strain with each BOX pattern was selected for genome sequencing to determine the Burkholderia species they belonged to. The results suggested that many of the isolates were clonal and that an internal outbreak might have happened at HGMGG.

Whole-genome analysis
The genomes of 12 selected strains were analyzed with the ANI and dDDH tools, considering that cutoff values of 95-96 and 70%, respectively, as the delineation of a species [40].  pattern XI as B. vietnamiensis (ANI 98.46% and dDDH 90.50%). The strain with the BOX pattern XII (strain 500H) was close to B. cenocepacia, but, according to the genome comparative analysis, it belongs to "B. servocepacia" (ANI 97.88% and dDDH 83.00%). However, this bacterial species has not yet been validated, and therefore it is not present in the LPSN. The genome sizes of the clinical strains were in the range of 6.5-8.8 Kbp, which is a genome size expected for bacteria in the Burkholderia genus. More genome features are presented in S2 Table. A phylogenomic analysis was performed with concatenating 92 genes with the UBCG pipeline. In this analysis, the type strain and four to five reference strains from each of B. cepacia, B. cenocepacia, B. multivorans, B. contaminans, B. vietnamiensis, B. pseudomallei and "B. servocepacia" species were included. The reference strains were selected on the basis of sharing more than 70% similarity by dDDH to the type strain to avoid using incorrectly classified species. The topology of the tree showed that all strains analyzed in this study were grouped with the respective species (Fig 2) and the strain 500H was clustered with "B. servocepacia" strains.
A dendrogram based on maximum parsimony from analysis of core SNPs for strain 294H was carried out comparing with a number of B. pseudomallei strains that were isolated from patients with melioidosis living in different locations [41]. Since strain 294H had the same BOX pattern with the rest of the B. pseudomallei strains analyzed in this study, this strain was randomly selected for genome sequencing and for SNPs analysis. The result showed that strain 294H clustered with B. pseudomallei strains that were isolated from people living in the US who had a history of travelling to Mexico and other countries in the American continent (Fig  3). Strain 294H was close to a strain isolated from a person with melioidosis living in Illinois, US, who previously had travelled to Mexico. This result shows the presence of the etiologic agent of melioidosis.

MLST analysis
The strains' genomes were analyzed in the PubMLST Database to determine whether they correspond to a defined sequence type (ST) or they represented a new ST. The B. cepacia strains belonged to a single ST, which was ST 9 (S2 Table).  Table). Since Burkholderia sp. 500H represents a novel species, the ST will be established once the species is proposed. These results showed that new STs were identified in some of the strains in this work.

Antibiotic susceptibility
Analysis with the VITEK2 system showed that all strains were resistant to ticarcillin/clavulanic acid, 87% were sensitive to trimethoprim/sulfamethoxazole, 94% were sensitive to meropenem, 82% were sensitive to levofloxacin and only 28% were sensitive to ceftazidime. These results showed that most of the strains analysed here are susceptible to elective antibiotics, although a small proportion are resistant, underscoring the relevance of defining the antibiotic sensitivity of all Burkholderia isolated in clinical settings.

Discussion
Bcc and B. pseudomallei are widely distributed in the world, but the evidence of their presence in Mexico is scarce, both in clinical settings and in the environment. Identification of these microorganisms at the species level is very problematic. Species from the Bcc can share up to 100% similarity in 16S rRNA gene sequences [39], prompting to the use of more efficient tools, such as MLSA and ultimately genome sequencing. Our analysis of Burkholderia strains isolated from human infections in two different hospitals in Mexico City showed that the identification system is poor, since all the strains were identified as B. cepacia, but we found that many belong to other Burkholderia species and even other genera. The proper identification of a bacterial species is important to provide the appropriate treatment, which is relevant in Burkholderia species as they are highly resistant to many antibiotics. For example, all strains studied were resistant to ticarcillin/clavulanic acid, but most of the strains were sensitive to meropenem. B. pseudomallei is efficiently treated with ceftazidime, meropenem or trimethoprim/sulfamethoxazole, and the B. pseudomallei strains studied here were sensitive to these antibiotics as well. Accurate identification of B. pseudomallei might lead to successful treatment of the patient. Unfortunately, this this did not happen in two fatal cases in northern Mexico; since the melioidosis caused by B. pseudomallei is barely known in Mexico, correct antibiotic treatment was not provided to the children, who died within a few hours after the admission to a hospital [20].
The analysis of 16S rRNA gene sequence showed poor resolution for the position of the Burkholderia strains in the phylogenetic tree. This is expected in this genus and with this gene. But when the atpD gene sequence was added to the analysis, the resolution of the bacterial positioning improved significantly, and the strains from HGMGG and HIMFG formed different clusters. Only the strains from HGMGG clustered with B. cepacia. This species is not typically recognized with highly distribution in clinical settings, as B. cenocepacia and B. multivorans are [5].
To group the strains and select some of them for genome sequencing, BOX-PCR was used, revealing 12 BOX patterns. One strain from each profile was selected for genome sequencing. The ANI and dDDH results confirmed that the strains from HGMGG belonged to B. cepacia, while the strains from HIMFG were identified as B. contaminans, B. multivorans, B. vietnamiensis and B. pseudomallei. A novel species in the Bcc was also found among the strains from HIMFG. In Mexico, B. cepacia has been found in several hospitals; however, the identification methodology is not always mentioned, and when automated systems, as used in this study, might not accurately identify the strains [7,11]. Although most of the strains from HGMGG were confirmed as B. cepacia, other strains had been mistakenly identified as this species, but belonged to other genera, such as Pseudomonas, Stenotrophomonas, and Providencia. B. cenocepacia and B. multivorans are Bcc species largely distributed in clinical settings, mainly affecting immunosuppressed patients [42,43,44,45,46,47]. In Mexico, B. multivorans strains have not been reported in clinical settings; this is the first report describing this species in Mexico. B. contaminans is considered an emerging CF pathogen in several countries [6], but not in Mexico. B. vietnamiensis strains have been isolated previously in Mexico, but from the environment. A number of these strains were obtained from maize, coffee and sorghum rhizosphere and rhizoplane in different geographic regions of Mexico [48]. B. vietnamiensis has been proposed as a plant growth-promoting rhizobacterium (PGPR) [49], but its use in agriculture is ethically incorrect since it is an opportunistic pathogen of the Bcc and has been found in patients with CF [50]. The presence of B. vietnamiensis in Mexican soils may represent a reservoir for the infection with this microorganism in immunosuppressed patients. Scientific and agricultural communities around the world are encouraged to stop proposing the use of Bcc strains as a PGPRs, as they are opportunistic pathogens.
By using genome sequencing, we corroborated the identification of strain 294H as B. pseudomallei [12]. This species is not well known in Mexico as few cases have been detected since 1958 [17]. The disease produced by B. pseudomallei, melioidosis, is present in Mexico, but it is a neglected disease because it is not accurately identified and reported. Recent cases with fatal outcomes [19,20], suggest that greater awareness of this bacterium and the disease is needed to prepare the health care staff to deal with the disease. Since the distribution of B. pseudomallei is worldwide, an increasing interest has emerged in the origin of B. pseudomallei strains. Several methods have been used to explore this origin, such as MLST and SNPs [38]. Some studies using SNPs have shown the association of clinical B. pseudomallei strains with environmental isolates, establishing an epidemiological link to understand the source of infection [51,52,53]. Our analysis with SNPs showed that strain 294H grouped with a strain that was isolated from a melioidosis patient in the US (Illinois) who had travelled to Mexico. Originally, the strain IL2014 formed a group with strains CA2009, MX2013, MX2014 and TX2015, all associated with clinical cases in patients residing in or visited Mexico [41]. However, strain 294H split this group and clustered closely with IL2014. This result show that more studies about the isolation of more B. pseudomallei strains are necessary to understand the origin of melioidosis infections both in Mexico but also in the US, since many US citizens travel for pleasure or work to Mexico. Moreover, since the US seems not endemic of melioidosis and many strains were isolated from people travelling to Mexico [41], our country might be endemic for B. pseudomallei.
All Burkholderia genomes were analyzed in the MLST database to determine whether they belonged to known STs or represented new ones. This database includes more than 100 different microbial species and genera, and it is intended to address several functional questions, such as activities of different variants that lead to key phenotypes. The Bcc PubMLST database functions to study the epidemiology and diversity of Bcc species. The B. cepacia strains from this study belonged to a single type, which was ST 9. This ST does not seem to be widely distributed since the database contains only two strains identified as ST 9: one isolated in Canada in 1997 (abdomen from non-CF patient) and a second isolated in Spain in 2009 (sputum from a CF-patient). The strains identified as B. contaminans were placed within ST 102 and ST 482. These STs are better represented, including 39 strains and 23 strains, respectively. The strains from B. multivorans and B. pseudomallei represented new STs; the former was assigned ST 1867 and the latter as ST 1872. However, the number of STs might be underestimated since not all Bcc species are present in the PubMLST database and the distribution of Bcc in the environment and clinical setting is poorly known in Mexico. However, we believe that the Mexican ST's might be endemic of Mexico, according to our previous analysis [54], but more studies should be carried out.
Additionally, two strains close to B. cenocepacia were found. By analyzing the genome of strain 500H, which shared the same BOX pattern with strain 501H, we found that it belonged to the novel species "B. servocepacia" [2]. The proposal of "B. servocepacia" was based on a comparison among many B. cenocepacia strains. The result of this analysis showed that B. cenocepacia could be divided in two major groups, and a few other strains might still represent novel species. Therefore, "B. servocepacia" was proposed and strain TAtl-371 was selected as the type strains. However, the species name of "B. servocepacia" has not been validated following the rules of the International Code of Nomenclature of Prokaryotes [55] and is not listed at the LPSN. Thus, we decided to maintain the name Burkholderia sp.
According to the identification of the strains, and taking into account the hospital of origin, the isolation source of the strains and the age of the person from whom the strains were isolated, it seems that no correlation exists. Only B. cepacia was identified from HGMGG but different species were obtained from HIMFG. The range of age was different in the patients from each hospital, but the isolation source might be important since in HGMGG most of the strains were isolated from bronchus secretions and in HIMFG from pharyngeal exudates.
Considering these results, the potential diversity of Burkholderia species in Mexico is not yet well understood. The distribution of species from the Bcc in Mexico is barely known, both in clinical settings and in the environment. The latter is important as a reservoir for potential infection in immunosuppressed patients which could explain the presence of Bcc and B. pseudomallei in human infections. The correct identification of Burkholderia species might help to provide accurate treatment for immunosuppressed patients.