Herbicide contaminated soil for enrichment studies was collected from Bien Hoa airbase (10°58′14.3″ N 106°48′19.3″ E), Dong Nai Province, Vietnam (Nguyen et al., 2021). Ten grams of herbicide contaminated soil was added in basal salt medium (BS, 40 ml per culture) containing KH₂PO₄ 0.5 g/l, (NH₄)₂SO₄ 0.25 g/l, MgSO₄ 0.2 g/l, CaCl₂ 0.5 g/l, and NaNO₃ 0.4 g/l (van der Zaan et al., 2012), and supplemented with 100 mg/l 2,4-D and 100 mg/l 2,4,5-T. Both compounds, more than 95% pure, were purchased from Sigma. Cultures were incubated in the dark on a rotary shaker at 30°C and 200 rpm (the first enrichment). After cultivation for 26 days, these cultures were then used to inoculate a second enrichment with 200 mg/l 2,4-D and 100 mg/l 2,4,5-T and further incubated for 25 days. Colony-forming units (CFUs) were isolated on nutrient agar (NA) at 30°C by serial dilution of the cultures at different time points after the start of the experiment. Nutrient agar medium contained 1.0 g peptone, 1.0 g meat extract, 0.5 g NaCl, and 1.5 g agar per 100 ml. Bacterial strains that were isolated as dominant CFUs were pre-cultured overnight in nutrient broth (NB) medium. The cells were harvested by centrifugation at 10,000 rpm for 5 min, washed two times with sterilized BS, and transferred into BS medium supplemented with 200 mg/l 2,4-D and 100 mg/l 2,4,5-T as inoculum for the degradation experiment under the same conditions as those for the enrichment cultures. A flask with dead cells, which were killed by heating the flask at 121°C during 15 min, amended with the herbicides was used as a control to display the effects of abiotic chemistry or denatured proteins. All experiments were independently performed in triplicate. Cell growth was determined by measuring the absorbance of cultures at 600 nm using a spectrophotometer (Multiskan GO, Thermo Scientific, U.S). We realize that the use of optical density values alone to detect growth of the bacterial cultures is less precise than other methods such as direct measurement of bacterial cell dry weight (DW) or flow cytometry. The concentrations of 2,4-D and 2,4,5-T in all culture samples were determined in time using tandem liquid chromatography-mass spectrometry (LC-MS/MS) (Aligent Technologies, 1200 Series). All experiments were performed according to the Biosafety Manual from the Vrije Universiteit Amsterdam entitled “Handling hazardous biological agents, genetically modified organisms and quarantine organisms” and approved by the biosafety officer of the same university.

The extraction and determination of the concentrations of 2,4-D and 2,4,5-T by LC-MS/MS along with the data analyses were performed exactly as described earlier (Nguyen et al., 2021).

The bacterial isolates were identified based on 16S rRNA gene sequence analysis. DNA of bacteria was extracted using the MoBio PowerSoil® DNA Isolation kit (Carlsbad, CA, USA) according to the manufacturer's instructions. Amplification was carried out with universal primers 8F (5′-AGAGTTTGATYMTGGCTCAG-3′) and 1512R (5′-ACGGYTACCTTGTTACGACTT-3′) as described previously (Weisburg et al., 1991). Reactions were performed in a Thermocycler (Biometra, Analytik Jena, Germany). PCR products were purified and sequenced by Macrogen Europe. These sequences were compared with known 16S rRNA gene sequences deposited in the GenBank database using the BLAST search at the National Center for Biotechnology Information (https://blast.ncbi.nlm.nih.gov/Blast.cgi).

Genomic DNA was extracted using the MoBio PowerSoil® DNA Isolation kit (Carlsbad, CA, USA). The quality and quantity of the extracted genomic DNA were checked by agarose gel electrophoresis and Qubit dsDNA HS Assay kit (Thermo Fisher Scientific, cat. no. Q32851). Genome sequencing was carried out by BaseClear B.V (The Netherlands) according to their in-house protocol. Single-end or paired-end sequence reads were generated using the Illumina NovaSeq 6000 or MiSeq system. The sequences generated with the MiSeq system were performed under accreditation according to the scope of BaseClear b.v. (L457; NEN-EN-ISO/IEC 17025). When paired-end sequencing is being performed, the “Number of reads” noted in Supplementary Table 2 is referring to read pairs. FASTQ read sequence files were generated using bcl2fastq2 version 2.18. Initial quality assessment was based on data passing the Illumina Chastity filtering. Subsequently, reads containing PhiX control signal were removed using an in-house filtering protocol, and reads containing (partial) adapters were clipped (up to a minimum read length of 50 bp). The second quality assessment was based on the remaining reads using the FASTQC quality control tool version 0.11.5. Supplementary Table 3 shows essential data of the sequencing and assembly of the genomes of the four selected isolates.

The raw Illumina reads were used directly for genome assembly using SPAdes (version 3.13.0), which is a De Bruijn graph assembler for short reads with a read error correction (Bankevich et al., 2012). The following SPAdes settings were used: –careful, -kmer: 55,77,127. Assembly quality was determined using QUAST (Gurevich et al., 2013) where we determined the number of scaffolds, total genome length, and N50 as a measurement for genome completeness. We found that the removal of scaffolds smaller than 5,000 bps did not significantly impact the total assembled genome length. Consequently, these scaffolds were removed in order to better facilitate the annotation. The scaffold size-adjusted genomes were annotated using several programs and databases. Gene prediction and initial annotation were done using Prokka (Seemann, 2014). Predicted genes were further annotated using Interproscan (Pfam and hamap databases) and the CAZy carbohydrate activated enzyme database (all families) (Jones et al., 2014; Lombard et al., 2014). Signal peptides were predicted with SignalP and secondary metabolite clusters using AntiSMASH (Weber et al., 2015; Almagro Armenteros et al., 2019). The closest homologs were determined by BLAST analysis using the Clusters of Orthologs Groups (COG) database for each predicted gene separately and by taking the species with the most hits as closest species (Tatusov et al., 2000). The NCBI-annotated genomes of BT1 2.2 (PRJNA586721), BT1 5.2 (PRJNA586745), BT1 9.2 (PRJNA586752), and BT1 10.2 (PRJNA587043) were deposited in GenBank with their BioProject numbers listed in between brackets.

The amino acid sequences of TftA homologs from the four isolated bacteria were compared with reference sequences available in the NCBI and KEGG databases. The sequences were aligned and a phylogenetic tree (Figure 3) was visualized using the MEGA version 7 software by the neighbor-joining method with 1,000 bootstrap replicates (Kumar et al., 2016).

The NCBI-annotated genomes of BT1 2.2 (PRJNA586721), BT1 5.2 (PRJNA586745), BT1 9.2 (PRJNA586752), and BT1 10.2 (PRJNA587043) that support the findings of this study were deposited in GenBank with their BioProject accession numbers listed in between brackets. Links to these data, initially for review only, are listed below. After acceptance of the manuscript, the data will become publicly available along with an assigned DOI number.

Pseudomonas aeruginosa strain: BT1 2.2: https://dataview.ncbi.nlm.nih.gov/object/PRJNA586721?reviewer=q7ko9k52t4a6uqlcrjdl3sdci2

Sphingomonas histidinilytica strain: BT1 5.2: https://dataview.ncbi.nlm.nih.gov/object/PRJNA586745?reviewer=i6f5dcie0oc5aamc0lu4h8g3ea

Bordetella petrii strain: BT 1 9.2: https://dataview.ncbi.nlm.nih.gov/object/PRJNA586752?reviewer=anls0i9ag2l5ehv9f3q7clktoi

Achromobacter xylosoxidans strain: BT1 10.2: https://dataview.ncbi.nlm.nih.gov/object/PRJNA587043?reviewer=ca8ottdc5qkfr9gdon77qj5ffi

Fourteen different types of bacterial colony were isolated from herbicide contaminated soil communities enriched in cultures with 2,4-D and 2,4,5-T as the sole sources of carbon and energy. The colonies of these strains on NA plates varied in shape and color. The results of 16S rRNA gene sequencing indicated that the dominant isolates were from genera Rhodanobacter, Arthrobacter, Pseudomonas, Achromobacter, Sphingomonas, Bordetella, Bosea, Pedobacter, Sphingopyxis, Staphylococcus, and Ancylobacter (Table 1). Among these isolates, four Gram-negative strains that showed up as the most dominant species in the enrichments, Pseudomonas, Achromobacter, Sphingomonas, and Bordetella, were selected to study their herbicide biodegradation capacities and their corresponding genetic make-up. BLAST results showed that the 16S rRNA gene sequences of the isolates BT1 2.2, BT1 5.2, BT1 9.2, and BT1 10.2 were highly similar (>99% identity) to those of P. aeruginosa, S. histidinilytica., B. petrii, and A. xylosoxidans, respectively. Therefore, we classified them as P. aeruginosa BT1 2.2, S. histidinilytica BT1 5.2, B. petrii BT1 9.2, and A. xylosoxidans BT1 10.2, respectively.

We set out a growth experiment where the isolates were cultured in basal salt medium supplemented with 2,4-D and 2,4,5-T as the sole sources of carbon and energy. We determined their growth curves (Figure 1A) and the removal of 2,4-D and 2,4,5-T in time (Figures 1B,C). Pseudomonas aeruginosa BT1 2.2 showed poor growth with a small increase in biomass after 7 days of cultivation, which subsequently decreased again after 10 days, most likely as a result of cell death and lysis. It did, however, show complete removal of 2,4-D and 2,4,5-T within 10 days of cultivation, but we have no data on the products that were formed. The other three strains showed substantial growth and complete removal of the herbicides within 10 days of cultivation. The fastest growers were B. petrii BT1 9.2 and A. xylosoxidans BT1 10.2, but the latter one had a lag time of 2 days before it started growing. Also, its final biomass was slightly less than that of B. petrii BT1 9.2. In both cases, maximum growth was reached after 4 days of cultivation. Sphingomonas histidinilytica BT1 5.2 needed a bit longer to reach its maximum growth yield. Due to its apparent slower growth rate, maximum biomass was obtained after 7 days, but the final yield was comparable to that of B. petrii BT1 9.2. The latter was the best performer concerning rates and yields of growth along with its ability to completely degrade 2,4-D and 2,4,5-T. To get a more comprehensive view of its capacities, we studied the responses of B. petrii BT1 9.2 in cultures with only 2,4-D or only 2,4,5-T and compared the growth data (Figure 2A) and removal of the herbicides (Figures 2B,C) with those from the culture with a mixture of both herbicides. The fastest growth and highest yields were observed in the cultures with only 2,4-D as the sole source of carbon and energy. Both growth rates and yields were slightly less in the mixture of the two compounds, and the poorest performance concerning rates and yield was observed with 2,4,5-T only. These growth properties coincided in part with the degradation rates of 2,4-D and 2,4,5-T. The former compound was removed within 2 days in the cultures with only 2,4-D, while it took 10 days in those with the mixture. Degradation of 2,4,5-T took 4 days of cultivation, while in the mixture this required 10 days.

We were also interested to understand the stoichiometry of conversion of the herbicides into biomass. For that, we made use of a few relevant numbers related to Escherichia coli taken from the Bionumbers archive (Milo et al., 2010) and assumed similar values for our isolated strains. At an OD₆₀₀ of 1, the number of E. coli cells is around 10¹²/L (Bionumbers ID 100985) with an average dry weight of 0.36 g/L (Bionumbers ID 109837). This is in the same order as when we approach these numbers from the average DW of an E. coli cell, which is 3.10⁻¹³ g DW (Bionumbers ID 100008), and, hence, 0.3 g/L DW in a culture with 10¹² cells. For the actual ratio of biomass formation over the catabolic branch one has to consider the percentages carbon in 2,4-D and in that of the biomass, which are 43 and 49% (assuming a biomass formula C₁H_1.8O_0.5N_0.24; Bionumbers ID 101800), respectively. Using the value of 0.36 g DW/L for a culture with an OD₆₀₀ of 1, we applied the calculations below. The increases in OD₆₀₀ of B. petrii BT1 9.2, S. histidinilytica BT1 5.2, and A. xylosoxidans BT1 10.2 after degradation of a total of 300 mg/L herbicides 2,4-D and 2,4,5-T are 0.33, 0.29, and 0.23, respectively. These values correspond to 119, 104, and 83 mg DW/L, from which one could calculate the formation of 0.40, 0.35, and 0.28 mg DW/mg 2,4-D mineralized, respectively. When we consider the carbon content in the herbicides and the DW, these numbers are 0.46, 0.40, and 0.33 mg C in DW/ mg C in 2,4-D, respectively.

Our data are also in compliance with data from the literature. The group of Tiedje grew an isolated Bradyrhizobium strain HW13 in media with 2,4-D as the sole source of carbon and energy and the growth curve showed an increase in the optical density of around 0.09 (32.4 mg DW/L) upon full mineralization of 125 mg/L 2,4-D (Kamagata et al., 1997). This corresponds to 0.26 mg DW/mg 2,4-D mineralized. The value calculated from the growth curve of C. pinatubonensis JMP134 growing on 3-chlorobenzoate (Pérez-Pantoja et al., 2000; Plumeier et al., 2002) is 0.46 mg DW/mg 2,4-D mineralized. Both values are in the same range as those of our isolates.

We inspected the four genomes for their genetic potential to degrade 2,4-D and 2,4,5-T to benzoate and further on to succinyl-CoA and acetyl-CoA. To this end, we gathered the amino acid sequences of all the enzymes identified in the KEGG database as being involved in these metabolic pathways. A list of our selected enzymes along with their properties and some additional parameters is shown in Supplementary Table 1. The KEGG database amino acid sequences of the selected genes were then used to look for similar sequences in the genome sequences of our four isolates. The script for these analyses is described in the materials and methods section. The results of these analyses for each organism are listed in Supplementary Table 3. Percentages identity between the KEGG protein sequences and those expressed from the selected genomes are listed as well. The threshold for selection of the orthologous proteins was set at a minimum of 45% identity. Importantly, all four species showed (initial) degradation of 2,4-D and 2,4,5-T so they should at least have a tfdA gene to degrade 2,4-D into 2,4-DCP, and a tftAB gene set to encode the dioxygenase for converting 2,4,5-T into 2,4,5-TCP, or as yet unknown enzymes that could take over that role. A phylogenetic tree of TftA homologs from the four isolates is shown in Figure 3 along with the reference enzymes from the KEGG and NCBI databases. These reference enzymes also include CadA, BenA, and AntA as we had reason to believe that some of the homologs that we found group closer to these proteins than to the TftA group. As such we were able to judge the minimum potential of the four species to degrade 2,4-D and 2,4,5-T and their further metabolites ultimately to succinyl-CoA and acetyl-CoA. We first constructed the reference pathways for degradation of the herbicides using the KEGG database maps “chlorocyclohexane and chlorobenzene degradation” and “benzoate degradation.” A graphical representation of these pathways is shown in Figure 4. We then reconstructed the pathways in our four isolated species by identifying their sequences orthologous to those of the reference pathways (Supplementary Figure 1). Below is a further description of the four isolates with emphasis on their genetic potential to convert 2,4-D and 2,4,5-T.

Achromobacter xylosoxidans BT1 10.2 has its tfdCDEF genes clustered along with a regulatory tfdS gene located upstream and divergently transcribed from it (Figure 5A). Downstream are genes encoding the subunits of a ring hydroxylating dioxygenase most similar to a 2-aminobenzenesulfonate 2,3-dioxygenase first described in an Alcaligenes species (Mampel et al., 1999). The catalytic part is composed of AbsAa and AbsAb and the electron donor is the flavodoxin, AbsAc. This gene organization and corresponding amino acid sequences are almost identical to those of B. petrii DSM 12804. This high degree of identity may be the result of a recent gene transfer event between ancestors of these species. Its tfdB gene is located elsewhere on the genome (ORF2511; Supplementary Table 3). We did not find an obvious tfdA gene nor tftAB genes for the initial conversion of 2,4-D and 2,4,5-T, respectively. A likely candidate for tfdA is ORF 3265 (Supplementary Table 3, 35% identity with the KEGG reference protein from C. pinatubonensis). Those for tftAB are ORFs 2313 and 2314 (around 32% identities with their KEGG counterparts from B. cepacia). Other genes potentially involved in further degradation of 2,4,5-T breakdown products are hbqR and chqB expressing HbqR and ChqB for the conversion of 2-HBQ via BT to 2-MA. Yet, the genes encoding the enzymes to convert 2,4,5-T into HBQ were not identified in this species.

Bordetella petrii BT1 9.2 has a full complement of tfdABCDEF genes encoding the enzymes for the sequential conversion of 2,4-D into 2-MA (Figure 5A). They are clustered and transcribed in the order tfdCDEFB with a lysR-type tfdS homolog upstream and divergently transcribed from that. A bit more downstream separated by four genes and with the transcription direction opposed to the other structural genes is the tfdA gene. One of the four genes sandwiched by the tfdA and remaining tfd genes is a gene encoding a PcaF-type enzyme responsible for cleavage of 3-oxoadipate-CoA into succinyl- and acetyl-CoA. The identity of the tfdBCDEF encoded proteins with the KEGG reference proteins from C. pinatubonensis JMP134 (accession no. M35097.1) is remarkably high, all of which ranging from 99 to 100%. Also, the gene organizations are identical suggestive of some kind of horizontal gene transfer between ancestors of these two Betaproteobacteria during recent evolution. We also compared the genetic organization of the isolated Vietnamese species with the one with its genome sequence deposited in the databases, B. petrii DSM 12804 (accession no. AM902716.1). Their gene organizations are different (Figure 5A). Also, the pairwise comparisons show that the identities are much lower than those between B. petrii BT1 9.2 and C. pinatubonensis (Supplementary Table 4). A second copy of the tfdA gene, not clustered with other tfd genes and with the expressed protein having 45% identity to the reference TfdA, was found elsewhere on the genome (ORF 3407). It has also the hbqR gene so it may be able to convert 2-HBQ to BT. However, it lacks reasonable homologs of enzymes to produce the substrate or to consume the product. Also, the hbqR gene is not clustered with other genes involved in 2,4,5-T degradation We did not find genes with high identity to tftAB genes but as this strain has the potential to perform the initial ring hydroxylating dioxygenase step, we assume that proteins related to the reference TftAB proteins are responsible for this reaction. Likely candidates to take over the role of TftA, all with identities to the reference TftA close to 30% are ORFs 524, 1633, 2507, and 3859. Their genes are all accompanied by homologs of tftB.

Sphingomonas histidinilytica BT1 5.2 has the genetic potential to express all enzymes for degradation of 2,4-D to succinyl-CoA and acetyl-CoA. Importantly, it does not have an apparent tfdA gene. The function of the tfdA gene product is likely taken over by a ring hydroxylating dioxygenase encoded by a cadABCD gene cluster that is sandwiched between the other tfd genes (Figure 5B). A similar gene organization is seen in Sphingomonas sp. TFD44 (accession no. AY598949.1), Sphingobium herbicidovorans MH (accession no. AJ628861.1), and Sphingomonas sp. ERG5 (accession no. KF494257.1) (Figure 4C) (Nielsen et al., 2013, 2017). Pairwise comparison of the gene products of the latter with those from our isolated strain indicates that they share a high degree of identity (more than 95%) both on the level of gene cluster organization and on the level of their amino acid sequences (Supplementary Table 4). Also, both clusters contain pcaIJF genes to make the enzymes for the conversion of 2-MA, the product of the Tfd proteins, into succinyl-CoA and acetyl-CoA, hence the potential of making a full pathway from 2,4-D to TCA-cycle intermediates. In addition to the gene cluster with tfd and cad genes, it has another one with copies of tfdCDEF genes on it, along with another set of pcaIJF genes. This one may be switched on for growth on chlorocatechols. Importantly, the suggested tfdD and tfdE genes in both clusters share a low identity with their counterparts from the KEGG reference proteins (in between 32 and 42%, see Supplementary Table 3). Nevertheless, we regard them as tfd genes because of their occurrence in the gene clusters with the complementary tfd genes. Its versatility in degrading chlorinated aromatics is further exemplified by the observation that it has a catABCD gene cluster under the apparent control of a catR regulatory gene (Figure 5C). Its gene organization looks quite similar to the one of P. aeruginosa BT1 2.2. Pairwise identities of the CatRABC proteins from both species are around 50% (Supplementary Table 4). Upstream of the cat gene cluster are two sets of genes that resemble those encoding the two subunits of the family of ring hydroxylating dioxygenases like AntAB, BenAB, and TftAB. It is tempting to speculate that these enzyme systems convert their as yet unknown aromatic substrates into catechol for further degradation to maleylacetate by the Cat proteins. It does not have obvious tftAB genes, yet it can convert it at least to 2,4,5-TCP. Hence another enzyme must be encoded by S. histidinilytica BT1 5.2. It might be that the CadABCD proteins can do that as they belong to the same family of ring-hydroxylating dioxygenases. Alternatively, that reaction may be carried out by the gene products of ORFs 5533 and 5532, which have around 30% identity with their counterparts TftA and TftB, respectively, from B. cepacia that serve as the KEGG database reference proteins.

Pseudomonas aeruginosa BT1 2.2 has catAB genes potentially involved in the 2,4-D breakdown. Their genes are clustered. A map of this gene cluster is in Figure 5C. A closer look at the gene cluster encompassing the catAB genes shows some interesting features as they make part of a typical catABC gene cluster, the proteins of which are involved in catechol degradation to 3-OXA. Transcribed in the same direction and sandwiched by the catB and catA genes is a gene homologous to catC and upstream but divergently transcribed a lysR-type regulatory gene possibly equivalent to catR (McFall et al., 1998; Nojiri et al., 2002). A similar gene organization is found in Pseudomonas putida mt-2 (Jiménez et al., 2014). Gene clusters composed of catRABC genes but organized differently are also found in R. erythropolis CCM2595 (accession no. FM995530.1) and Burkholderia sp. TH2 (accession no. AB035483.1). More upstream are two gene clusters encoding the subunits of ring-hydroxylating oxygenases, both of which under the apparent control of AraC-type regulators, their genes flanking these two gene clusters. The one most close to the cat gene cluster has two structural genes with the highest similarity with antABC genes to make an anthranilate 1,2-dioxygenase, AntAB as the catalytic core, and AntC as the flavodoxin-type electron donor. The one more downstream is most likely a benABCD gene cluster encoding a benzoate dioxygenase (BenABC) and a dihydroxy cyclohexadiene carboxylate dehydrogenase (BenD). Burkholderia sp. TH2 has also a set of benABCD genes adjacent to its cat gene cluster (Figure 5C). The phylogenetic tree with TftA homologs including AntA and BenA proteins shows their position in the sub-groups, which supports their suggested functions. Both enzyme systems, AntABC and BenABCD, convert their substrates, anthranilate, and benzoate, respectively, to catechol. Altogether, these findings suggest that P. aeruginosa BT1 2.2 is specialized in the degradation of benzoate and related compounds as it has ben and cat gene clusters for the sequential conversion of benzoate via catechol to 3-oxoadipate. Genes encoding PcaIJ and PcaF are found elsewhere on the genome for a further breakdown of 3-oxoadipate to succinyl-CoA and acetyl-CoA. It can convert 2,4-D at least into 2,4-DCP, and we speculate that this is achieved by the product of ORF3804, which has 38% identity with the KEGG reference tfdA from species C. pinatubonensis JMP134. An alternative gene set for tftAB genes might be ORFs 1030 and 1031. The protein from ORF 1030 has 30% identity with the KEGG reference TftA from B. cepacia AC1100.