All chemicals were purchased from Sinopharm Chemical Reagent Co. Ltd., including tetracycline standard, tryptone, yeast extract, NaCl, glucose, chromatography grade methanol and acetonitrile, analytical grade sodium dihydrogen phosphate, oxalic acid, disodium ethylenediaminetetraacetate, and citric acid.

The soil sample was obtained from the sludge in the sediment pool of a pig factory using tetracycline antibiotics as feed additive for a long-term (4 a) in Changsha. The average temperature and pH value of the sampling location in this study was 20°C and about 6.8. The latitude and longitude of our sampling locations are 28.1932 and 112.6973, respectively.

Lysogeny broth (LB) medium contained the following (in g L⁻¹): tryptone, 10; yeast extract, 5; NaCl, 5; glucose, 10. The pH of LB medium was adjusted to 7.0–7.2 before autoclaving at 121°C for 20 min.

0.1 g soil sample was made into suspension with 10 mL sterile water. One milliliter supernatant was added to the LB medium (An additional volume of tetracycline stock solution was added to the culture medium so that the tetracycline concentration was 20 mg L⁻¹). The mixed solution was incubated at 30°C on a dark shaker at 150 rpm. After 3 days, 0.5 mL of the solution was inoculated into the LB medium containing 40 mg L⁻¹ tetracycline and the other conditions remained unchanged, and then cultured for 6 days. The above procedures were repeated until the final concentration of tetracycline reached 120 mg L⁻¹. After the enrichment, the solution was serially diluted from 10⁻¹ to 10⁻⁹, and then 0.5 mL of each diluted solution was added to petri dishes and daub uniformly. Six colonies with different morphology were obtained after multiple purifications, and then inoculated separately into LB medium for cultivation. Only the isolate with a high level of degradation efficiency of tetracycline was further analyzed. Finally, Pandoraea sp. XY-2 was identified through 16S rRNA sequence analysis as previously described (Mangwani et al., 2015).

Concentration of tetracycline residues were analyzed by high-performance liquid chromatography (HPLC). A C18 reversed-phase column (Agilent Technologies) was operated at 30°C, with a 0.5 mL min⁻¹ mobile phase consisting of 10 mmol L⁻¹ sodium dihydrogen phosphate solution and acetonitrile. The injection volume was 20 μL, and the column gradient elution was monitored by a UV detector at 360 nm. Five milliliter solution from each reactor was centrifuged at 8,000 rpm at 4°C for 15 min. Four milliliter supernatant was added to a centrifuge tube containing 4 mL Na₂EDTA-Mcllvaine buffer, then 1 mL n-hexane and 1 mL chloroform were added, and the mixture was shaked for 2 min and sonicated for 30 s. After the mixture was centrifuged at 8,000 rpm for 15 min, the supernatant was filtered through a 0.22 μm microporous membrane filter and preserved at 4°C for further testing with the previously described chromatographic conditions.

In order to obtain the optimal conditions for Pandoraea sp. XY-2 biodegrading tetracycline, we studied the growth of Pandoraea sp. XY-2 and the biodegradation percentage of tetracycline under various culture time (i.e., 1, 2, …, 8, 9 day), temperatures (i.e., 20, 25, 30, 35, 40°C) and initial pH values (i.e., 5.0, 6.0, 7.0, 8.0, 9.0). Only biodegradation was considered because hydrolysis of tetracycline could be ignored during the testing period (Supplementary Figure 1). In these experiments, the amount of each bottle of culture solution was 150 mL. Pandoraea sp. XY-2 culture solution was inoculated to LB medium containing 50 mg L⁻¹ tetracycline at 3% inoculum concentration, and then placed on a dark shaker at 150 rpm. During the culture time tests, the temperature and initial pH were 30° C and 7.0, respectively. During the different temperature tests, the initial pH was 7.0. During the pH tests, the temperature was 30°C. The final growth of Pandoraea sp. XY-2 and the biodegradation percentage of tetracycline in temperature and pH tests were measured only at the 6th day. All experiments were performed in triplicates. The growth of Pandoraea sp. XY-2 was measured by a hemocytometer and the concentration of tetracycline in the medium was measured using HPLC.

Genomic DNA of Pandoraea sp. XY-2 was extracted using the Genomic DNA Extraction Kit (Tian Gen biochemical technology Co. Ltd., Beijing, China).

Genomic DNA was sheared in a Covaris g-TUBE, and DNA fragments were then damage repaired, terminal repaired and purified using AMpure PB magnetic beads. SMRTbell libraries prepared with DNA Template Kit 2.0 and quantified by Qubit concentration. Agilent 2,100 was used to detect the size of the inserted fragments. Then, sequencing carried out using the PacBio RS II SMRT sequencing technology (Beijing Nuohe Zhiyuan Technology Co., Ltd., Beijing, China).

De novo genome assembly of PacBio reads for XY-2 (91,372 reads totaling 1,072,300,065 bp, mean length 11,735 bp, N50 length 16,735 bp) was conducted using the hierarchical genome assembly process (HGAP3 protocol and SMRT Link v.5.0.1 software) (Rhoads and Au, 2015). The mapping of assembled contigs onto long reads from the PacBio platform was performed by BlasR (Chaisson and Tesler, 2012). After mapping, long reads with high accuracy were assembled by the Celera Assembler v8.0 (Science New York, 2000). Then, the assembly resulted in one polished contig, which was a circular chromosome with no gaps.

Gene prediction carried out by GeneMarkS software (http://topaz.gatech.edu/). The tRNAscan-SE v1.3.1 and rRNAmmer v1.2 were used for identifying tRNA and rRNA, respectively (Schattner et al., 2005; Lagesen et al., 2007). Additionally, online Web-based tool Island Viewer v.3 was used for detecting genomic islands (Dhillon et al., 2015). To predict Interspersed Repeats, RepeatMasker was used (Chen, 2004). Tandem Repeats were identified for by using Tandem Repeat Finder v.4.09 (Benson, 1999). Finally, the gene functions were further investigated by BLASTP using the following five different function databases: the carbohydrate-active enzymes (CAZy) database (http://www.cazy.org/), the non-redundant (NR) protein database (http://www.ncbi.nlm.nih.gov/protein), the clusters of orthologous groups (COG) database (https://www.ncbi.nlm.nih.gov/COG/), the GO database (http://www.geneontology.org/), and the KEGG database (http://www.genome.jp/kegg/pathway.html).

Before ANI analysis, the 16S rRNA sequence of type strains of Pandoraea spp. were retrieved from the EzTaxon server (http://www.eztaxon.org/, Jongsik et al., 2007). Comparative sequence analysis of the 16S rRNA gene of Pandoraea sp. XY-2 agaist the type strains of Pandoraea spp. was implemented within the server. The phylogenetic tree of Pandoraea sp. XY-2 based on the 16SrRNA gene sequences was constructed using the NJ method in MEGA v7.0 with 1,000 bootstraps replications (Kumar et al., 2016).

Reference genomes for comparison purposes were retrieved from the GenBank database (http://www.ncbi.nlm.nih.gov/genbank/). Whole genome comparisons to determine the ANI between genome sequences were produced using JSpeciesws database (http://jspecies.ribohost.com/jspeciesws) that use Blast alignments to evaluate whole genome homologies. Cut-offs for species delineation were 95% ANI on 69% of conserved DNA according to Goris et al. (2007). In addition, we also calculate the dDDH values using Genome-to-Genome Distance Calculator (Alsaari et al., 2015). The 70% species cut-off of dDDH usually kept in taxonomic studies of bacteria (Konstantinidis and Tiedje, 2007).

The other 35 Pandoraea strains (Table 1) from NCBI database were selected to perform the phylogenetic analysis of Pandoraea sp. XY-2. Firstly, the single-copy core genes were identified by core-pan analysis between Pandoraea sp. XY-2 and the 35 Pandoraea strains. The number and total length of concatenated single copy genes were 1,903 and 27,797 characters. Then, the multiple sequences alignment of single-copy core genes was performed using MUSCLE software. Finally, phylogenetic tree was generated using the NJ method in Mega v7.0 with 1,000 bootstraps replications.

We used the Bacterial Pan Genome Analysis tool (BPGA) pipeline (Chaudhari et al., 2016) to identify orthologous groups among Pandoraea testing genomes and to extrapolate the pan-genome models of applying default parameters. The statistics information and details source of all testing strains used are listed in Table 1 and Supplementary Table S1. We defined the set of genes shared among all strains as their core genome, the set of genes shared with more than two but not all strains as their accessory genome, the set of genes not shared with other strains in each strain as the specific genes and the set of non-homologous genes with all testing genomes as the pan genome. Finally, functional categories were assigned to gene clusters by the KEGG annotation previously described.

The size of Pandoraea pan-genome was extrapolated implementing a power law regression function P_s = κn^γ using a built-in program of BPGA pipeline (Chaudhari et al., 2016), where P_s is the total number of non-orthologous gene families within its pan-genome, n is the number of tested strains, and both κ and γ are free parameters (Tettelin et al., 2008). The exponent γ < 0 suggests the pan-genome is “closed” considering the size of the pan-genome is reaching a constant as extra genomes adding in. Conversely, species is predicted to harbor an “open” pan-genome for 0 < γ < 1. Furthermore, the size of the core-genome was extrapolated fitting into an exponential decay function F_c = κ_cexp^(−n/τc) + Ω, with a built-in program of BPGA pipeline (Chaudhari et al., 2016), where F_c is the number of core gene families, while κ_c, τ_c, and Ω are free parameters (Tettelin et al., 2005).

For variation analysis, we use P. apista DSM16535^T as reference strain. The single nucleotide polymorphisms (SNPs) were identified by aligning scaffolds to reference using MUMmer 3.06 package (Kurtz et al., 2004). Reliable SNPs were obtained after using BLAST, TRF, and Repeatmask software to predict the repeat region of the reference sequence, and filter out the SNPs located in the repeat region. InDels (insertion and deletion variations) were detected through gap extension alignment between the whole genome de novo assembly and the reference using LASTZ software (Harris, 2007). Synteny analysis was performed using MUMmer v3.06 package.

The coding sequences that are resistant to antibiotic was predicted using the Resistance Gene Identifier (RGI) software in the Comprehensive Antibiotic Resistance Database (CARD) (Mcarthur and Wright, 2015). After identifying the tetracycline resistance genes, each resistance gene was compared on NCBI database to find the genes with high similarity. Then, the evolutionary history of genes for tetracycline resistance was inferred using MEGA v7.0 (Kumar et al., 2016) with the Neighbor-Joining method with 1,000 bootstraps.

The sequence data of the whole genome of Pandoraea sp. XY-2 have been deposited to NCBI data base, with the accession number of CP030849.

After enrichment, isolation, and purification, 6 strains named XY-1, XY-2, XY-3, XY-4, XY-5, and XY-6 were obtained. The colony morphology of strain XY-1 was round, white, opaque, while that of strain XY-3 was dry, white, uneven edges. The colony morphology of strain XY-4, XY-5, and XY-6 was cotton-like, snowflake and fluffy shape, respectively. All the strains were capable of biodegrade tetracycline, and the biodegrade percentage were between 36 and 65%. The strain XY-2 was selected for further analyses due to its highest ability to biodegrade tetracycline. The colony of strain XY-2 on LB solid medium was smooth and moist, and was light milky white with neat edges. It is a gram-negative bacterium and has a short rod shape under a scanning electron microscope (Supplementary Figure 2) and has a size of ~0.7 μm × (1.2–3) μm. After the 16S rDNA sequence analysis of strain XY-2, it was concluded that strain XY-2 belongs to genus Pandoraea in molecular phylogenetic classification and it was named Pandoraea sp. XY-2.

Figure 1 shows the effect of incubation time on the growth of Pandoraea sp. XY-2 and tetracycline biodegradation percentage. The growth of Pandoraea sp. XY-2 entered the logarithmic phase at the 3rd day, the biomass reached 2.1 × 10¹⁰ cfu/mL at the 6th day, and then the growth entered the stationary phase. The biodegradation percentage of tetracycline increased significantly during the first 6 days and it reached 73.8% at the 6th day. After that, there was almost no significant change with prolonged incubation.

Effect of incubation temperature on the growth of Pandoraea sp. XY-2 and tetracycline biodegradation ability is shown in Figure 2. The growth of Pandoraea sp. XY-2 was greatly affected by temperature. With the increasing temperature from 30 to 40°C, the biomass of Pandoraea sp. XY-2 decreased to 1.57 × 10¹⁰ cfu/mL and the biodegradation percentage of tetracycline reached the maximum value and then remained stable at the 6th day.

The proper pH value can affect not only the growth rate of microorganisms, but also the activity of enzymes. In this study, the growth of Pandoraea sp. XY-2 was significantly inhibited when the initial pH was 5 and 9 (Figure 3). The Pandoraea sp. XY-2 grew most vigorously and the biodegradation percentage of tetracycline reached the maximum value at pH 7. It was concluded that the growth of Pandoraea sp. XY-2 and tetracycline biodegradation percentage generally showed a trend of increasing first and decreasing then with the increasing initial pH.

Combined with the growth of Pandoraea sp. XY-2 and tetracycline biodegradation percentage, we selected a pH of 7, temperature of 30°C and culture period of 6 day as the optimum culture conditions for subsequent experiments.

After sequencing, a total of 91,372 reads were obtained. Genome characteristics of the newly sequenced genome of Pandoraea sp. XY-2 are shown in Table 2. The Pandoraea sp. XY-2 has a single circular chromosome and the estimated genome size was 5,056,206 bp with a GC content of 63.76% (Figure 4). No plasmid was found in the genome sequence of this bacterium. The number of predicted genes was 4,974 with an estimated total length of 4,382,283 bp which makes up 86.67% of the entire genome. Out of all the predicted genes, 3,489 (70%) could be functionally annotated in the NCBI non-redundant (nr) database. A total of 1,485 hypothetical proteins were predicted for XY-2. There are 77 RNA genes consisting of 65 tRNA and 12 rRNA (4 5S rRNA, 4 16S rRNA, and 4 23S rRNA) genes. Eight GIs were detected in the genome of the Pandoraea sp. XY-2. Approximately 5,470 bp of interspersed repeats were found in the genome of the Pandoraea sp. XY-2, which accounted for 0.1082% of its genome size, whereas 27,270 bp of tandem repeats were found, which accounted for 0.5393% of its genome size. Tandem repeats were the predominant type of repeat region, accounting for 83.3% of the repeat regions in the Pandoraea sp. XY-2 genome.

Secreted CAZymes are crucial for bacteria biological activity. CAZymes are divided into five classes: glycoside hydrolases (GHs), carbohydrate esterases (CEs), polysaccharide lyases (PLs), glycosyltransferases (GTs), and auxiliary activities (AAs). Based on the CAZy annotation pipeline, we identified 65 CAZymes, including 5 AAs, 16 GHs, 30 GTs, 3 CEs, and 11 carbohydrate binding modules (CBMs) in the Pandoraea sp. XY-2 genome (Supplementary Figure 3). In the present study, CAZyme annotation revealed that Pandoraea sp. XY-2 contains a total of 2 families with 5 AAs in its genome sequence. Additionally, AA family classification also revealed that the components of AAs were 3 AA3 family members [glucose-methanol-choline (GMC) oxidoreductase, alcohol oxidase, aryl-alcohol oxidase/glucose oxidase, cellobiose dehydrogenase, pyranose oxidase), and 2 AA6 family members (1,4-benzoquinone reductase).

Clusters of orthologous groups (COG) database annotation showed that a number of genes are related to the basic functions of bacterial cells, including amino acid transport and metabolism (E), transcription (K), and energy production and conversion (C) (Supplementary Figure 4). Particularly, the number of genes participated in carbohydrate transport and metabolism (G), secondary metabolites biosynthesis, transport and catabolism (Q) and intracellular trafficking, secretion, and vesicular transport (U) were 244 (5.4%), 166 (3.6%), and 79 (1.7%), respectively. The GO annotation revealed that cellular processes and metabolic processes are very active in biological process; cell and cell part play a dominant role in cellular component; binding and catalytic activity play important roles in molecular function (Supplementary Figure 5). Furthermore, KEGG annotation showed that a total of 115 genes participated in the xenobiotic biodegradation and metabolism (Supplementary Figure 6 and Supplementary Table S2). It is worth noting that there were 5 genes encoding monooxygenase, 9 genes encoding dioxygenase, nine genes encoding glutathione S-transferase (GSTs) and 92 genes encoding other types of enzymes. Previous reports showed that tetracycline can be more readily biodegradable by using the detoxifying enzyme GSTs (Park and Choung, 2007). Accordingly, we selected all the 9 GSTs-encoding genes in Pandoraea sp. XY-2 together with GSTs-encoding genes from close-related strains for comparative analysis to shed light on the variation in the ability to degrade tetracycline between Pandoraea spp. The phylogenetic analysis of nine GSTs-encoding genes derived from Pandoraea sp. XY-2 and other close-related strains revealed that the GSTs-encoding genes were present in both environmental and clinical isolates and the evolution of these genes had no major variation (Supplementary Figure 7). We also calculated the GSTs-encoding genes contents in Pandoraea spp. and the results showed that the contents of these genes did not change significantly between environmental and clinical isolates (Supplementary Table S1).

The phylogenetic tree of Pandoraea sp. XY-2 based on 16SrRNA gene sequence was shown in Supplementary Figure 8. Pandoraea sp. XY-2 was most closely related to the P. apista LMG16407^T (99.59 similarity). The 16S rRNA gene of Pandoraea sp. XY-2 having a similarity of 96.9–99.72% compared to the type strains of Pandoraea spp. (Supplementary Table S3). (Yarza et al., 2008) have reported that a new proposed genus would have ~6% divergence in 16S rRNA gene sequence from its closest genus. Therefore, Pandoraea sp. XY-2 was not a new genus. Stackebrandt and Ebers (2006) suggested that a 16S rRNA gene sequence similarity of 98.7–99% should become the boundary for delineation of prokaryotic species. Above this value, there is a need to prove by genome to genome comparisons whether distinct isolates belong to the same species. Accordingly, genomes comparisons of Pandoraea sp. XY-2 against all Pandoraea spp. were carry out for species circumscription.

ANI is considered to be the most relevant comparative parameter used for bacterial species delineation (Rosselló-Móra and Amann, 2015). Using ANI, the bacterial species delimitation borders can be set to about 94–96% identity, which would generally represent ~70% dDDH (Richter and Rossellómóra, 2009). Table 3 contain the matrix of nucleotide identities between whole genomes calculated under JSpeciesws database. The ANI% values of Pandoraea sp. XY-2 against all reference strains were < 85% (Table 3), which was a value far below the threshold of 94–96% identity that would serve as a boundary for species circumscription. Besides, the dDDH results showed that the dDDH% values of Pandoraea sp. XY-2 against all reference strains ranging from 14.9 to 52.2% (Supplementary Table S4).

A phylogenetic analysis based on the single-copy core genes was conducted to determine the relationship of Pandoraea sp. XY-2 with other 35 Pandoraea strains. Except for the sequencing data of Pandoraea sp. XY-2, other whole genomic data of 35 Pandoraea strains were downloaded from NCBI. The Pandoraea sp. XY-2 was clustered together with Pandoraea sp. 64–18 and Pandoraea sp. B-6 (Figure 5), which was indicative of intimate intragenus relationship among Pandoraea sp. XY-2, Pandoraea sp. 64–18 and Pandoraea sp. B-6. In addition, Pandoraea sp. XY-2 branches near the base of the Pandoraea sp. 64–18 and Pandoraea sp. B-6. We inferred that Pandoraea sp. XY-2 was likely arised before the Pandoraea sp. 64–18 and Pandoraea sp. B-6.

A pan-genome for the Pandoraea sp. XY-2 and 35 other Pandoraea strains was determined by PGAP pipeline. The total pan-genome for the 36 compared Pandoraea strains encompasses 186,995 putative protein-coding genes. Of these, 1,903 (1.02% of total pan-genome) were core conserved genes across all 36 strain genomes compared. Out of these 186,995 genes, 2,280 (1.22%) were represented in the accessory genomes of Pandoraea sp. XY-2. After comparing strain-specific genes, the Pandoraea sp. ISTKB had the highest amount of these (n = 899), while the isolate sequenced in this study had 414 specific genes. The number of specific genes in all 36 strains ranges from 2 to 899 (Figure 6). The KEGG annotation of specific genes of Pandoraea sp. XY-2 showed that a total of 10 genes involved in the environmental information processing, including one TetR family protein-encoding gene tetR (XY.2_GM004468) and one DHA1 family protein-encoding gene tetG (XY.2_GM004469) (Supplementary Figure 9).

Previous reports showed that mathematical extrapolation of pan-genome would be highly reliable provided that sufficient genomes (>5) are involved (Vernikos et al., 2015). The deduced power law regression function [P_s (n) = 3119.72 n^0.50637] revealed that the pan-genome of Pandoraea had a parameter (γ) of 0.50637 falling into the range 0 < γ < 1, which indicating that the pan-genome was open (Figure 7) (Tettelin et al., 2008). In addition, the deduced exponential regression [F_c (n) = 3900.47e^{−0.0178922n}] revealed that the extrapolated curve of core genome followed a steep slope, reaching a minimum of 1,903 gene families after the 36th genome was added (Figure 6).

The other 18 Pandoraea strains, which genome sequences level were complete, from NCBI database were selected to perform the phylogenetic analysis of Pandoraea sp. XY-2. The results showed that Pandoraea sp. XY-2 was clustered together with P. apista DSM16535^T (Supplementary Figure 10). Combined with that it has not yet been reported that tetracycline could be biodegraded by Pandoraea spp., so we selected Pandoraea sp. XY-2 and P. apista DSM16535^T to perform a genome-wide comparison. The genome-wide comparison identified all the SNPs, a total of 372,574 high-quality SNPs for Pandoraea sp. XY-2 were obtained (Supplementary Table S5). Most of the high-quality SNPs in the Pandoraea sp. XY-2 were located in CDS regions (333,149, 89.42%), with only 39,425 (10.58%) located in intergenic regions. A total of 67,476 non-synonymous SNPs and 264,974 synonymous SNPs were located within CDS regions, which resulted in a non-synonymous/synonymous (dN/dS) ratio of 0.255.

Gene gains and losses have been considered as one of the most significant contributors to functional changes (Nei and Rooney, 2005). We analyzed the insertion/deletion variations (InDels) between Pandoraea sp. XY-2 and P. apista DSM16535^T and the results are shown in Supplementary Table S6. We identified 581 InDels in the chromosome of Pandoraea sp. XY-2 and the proportion of base deletion (60.76%) is higher than that of base insertion (39.24%). Most of the InDels were located in intergenic regions (388, 66.78%). Furthermore, synteny comparison of genomes between Pandoraea sp. XY-2 and P. apista DSM16535^T showed that the two genomes have 550 collinear blocks covering 73.91 % of Pandoraea sp. XY-2 genome. On the other hand, some regions between Pandoraea sp. XY-2 and P. apista DSM16535^T displays inversion, translocation, and Tran+Inver, revealing different syntenial relationships to the reference sequence (Figure 8).

From CARD analysis, a total of 218 genes were identified as having a role in the resistance to various antibiotics, of which 9 were tetracycline resistance genes (Supplementary Table S7). The type and number of tetracycline resistance genes were 6 tetA(48), 1 tetB(P), 1 tetG, and 1 tetT (Supplementary Table S8). Resistance is primarily due to either energy-dependent efflux of tetracycline or protection of the ribosomes from the action of tetracycline. TetA(48) and tetG belong to the efflux pump gene, while tetB(P) and tetT belong to the ribosome protection gene (Roberts, 1996). Previous studies have shown that efflux pumps play a major role in both organic contaminants tolerance and bioremediation (Fernandes et al., 2003), and combined with KEGG annotation results of specific genes in Pandoraea sp. XY-2, so we selected tetA(48), tetG, and tetR genes for further origin and evolution analysis.

The deviant GC content is used as a detect method of HGT (Xie et al., 2014). The average GC contents of the tetA(48) genes were higher than those of the average of the entire genomes (63.33–65.57 vs. 62.5–64.7) except that the average GC content of the tetA(48) genes in Pandoraea sp. XY-2 was slightly lower than that of the average of the genome (63.72 vs. 63.76) (Figure 9). It is worth noting that the tetG gene and tetR gene were specific genes of Pandoraea sp. XY-2, whose GC content were significantly lower than that of the genome (57.74 vs. 63.76 and 59.17 vs. 63.76, respectively). These data showing variation of GC content between tetracycline resistance genes and the corresponding genomes, which indicated that these tetracycline resistance genes may be acquired in Pandoraea strains by HGT.

To gain insights into the origin of tetracycline resistance genes in Pandoraea, three NJ phylogenetic trees were constructed based on the tetA(48), tetG, and tetR genes encoding protein sequences (Figures 10–12). As shown in Figure 10, phylogenetic analysis revealed that three tetA(48) genes (XY.2_GM000586, XY.2_GM000605 and XY.2_GM004567) of Pandoraea and tetA(48) genes of Paraburkholderia spp., Burkholderia spp., and Caballeronia spp. were sister groups. It is obviously to note that, the other three tetA(48) genes (XY.2_GM001703, XY.2_GM003133 and XY.2_GM003479) of Pandoraea sp. XY-2 formed monophyletic clade, respectively. Furthermore, the tetG gene from Pandoraea sp. XY-2 grouped to tetG gene of S. enterica, V. cholerae, P. aeruginosa, and Proteobacteria (Figure 11). The tetR gene from Pandoraea sp. XY-2 grouped to tetR gene of S. enterica, V. cholerae, A. baumannii, and F. marinus (Figure 12). These results imply that three tetA(48) genes (XY.2_GM000586, XY.2_GM004567, and XY.2_GM000605) may be acquired via HGT from members of Paraburkholderia, Burkholderia, and Caballeronia. While, the other three tetA(48) genes (XY.2_GM001703, XY.2_GM003133, and XY.2_GM003479) may be acquired via HGT from unidentified sources in early evolutionary history. The tetG gene may be acquired via HGT from S. enterica, V. cholerae, P. aeruginosa and Proteobacteria, and tetR from S. enterica, V. cholerae, A. baumannii, and F. marinus.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.