Subfunctionalization influences the expansion of bacterial multidrug antibiotic resistance
BMC Genomics volume 18, Article number: 834 (2017)
Antibiotic resistance is a major problem for human health. Multidrug resistance efflux pumps, especially those of the Resistance-Nodulation-Cell Division (RND) family, are major contributors to high-level antibiotic resistance in Gram-negative bacteria. Most bacterial genomes contain several copies of the different classes of multidrug resistance efflux pumps. Gene duplication and gain of function by the duplicate copies of multidrug resistance efflux pump genes plays a key role in the expansion and diversification of drug-resistance mechanisms.
We used two members of the Burkholderia RND superfamily as models to understand how duplication events affect the antibiotic resistance of these strains. First, we analyzed the conservation and distribution of these two RND systems and their regulators across the Burkholderia genus. Through genetic manipulations, we identified both the exact substrate range of these transporters and their eventual interchangeability. We also performed a directed evolution experiment, combined with next generation sequencing, to evaluate the role of antibiotics in the activation of the expression of these systems. Together, our results indicate that the first step to diversify the functions of these pumps arises from changes in their regulation (subfunctionalization) instead of functional mutations. Further, these pumps could rewire their regulation to respond to antibiotics, thus maintaining high genomic plasticity.
Studying the regulatory network that controls the expression of the RND pumps will help understand and eventually control the development and expansion of drug resistance.
Antibiotic resistance is a major challenge for the twenty-first century . A large number of bacteria resistant to multiple classes of antibiotics has emerged worldwide, making treatment difficult [2, 3]. Multidrug antibiotic resistance (MDR) is a highly attractive model to study the evolution of gene function since hypermutation, complex interrelationships between drug resistance and fitness, compensatory evolution, and epistasis affect how resistance evolves and spread . Although different mechanisms influence the emergence of a MDR phenotype, high-level intrinsic MDR, particularly in Gram-negative bacteria, stems from the combined action of efflux pumps and adaptive modifications of the cell envelope .
The structural components of MDR efflux pumps (MDR EP), usually chromosomally encoded, are evolutionarily conserved within species [reviewed in ]. Indeed, MDR EP are ancient elements in bacterial genomes, suggesting that their functions predate the resistance to antibiotics during the treatment of human infections . Bacterial genomes often contain multiple copies of different MDR EP classes, which have arisen by gene duplications [6, 7]. The importance of gene duplication in driving expansion and diversification of drug-resistance mechanisms and its role as an adaptive response to antibiotic treatment are well-known [8,9,10]. In general, extra gene copies after duplication events are redundant and free from selection pressure, as they do not add anything to the organism’s capacity to perform the original (duplicated) function . The fate of paralogous copies is strongly debated and no single model can include all the possible alternative scenarios . Duplicated genes rarely exhibit de novo functions (neofunctionalization); more commonly, the functions of the original gene(s) are split into multiple functions among the duplicate genes (subfunctionalization) . The split into subfunctions among the different paralogs may occur when the functions of the original gene(s), which were previously under unified genetic control, acquire their own independent regulation .
Among MDR efflux pumps, the Resistance-Nodulation-Cell Division (RND) superfamily is particularly intriguing . Members of this superfamily of proton/drug antiporters located in the inner membrane have several roles including bacterial virulence, quorum sensing, plant–bacteria interactions, and detoxification of metabolic intermediates and toxic compounds, such as heavy metals, solvents, or antimicrobials . The regulation of RNDs depends on global and/or local regulators, resulting in a multilayered control of gene expression in response to specific stimuli. The high efficiency of RNDs in extruding substrates is due to their associations with outer membrane channel (OMP) and periplasmic membrane fusion (MFP) proteins. These tripartite complexes extrude molecules to the extracellular milieu preventing their accumulation in the periplasmic space . The genes encoding the three protein components (RND, OMP, and MFP) are often organized in an operon and several different RND operons may be embedded in the same genome [3, 17]. These efflux pumps have been mostly studied in Gram-negative bacteria, such as Pseudomonas aeruginosa, Escherichia coli, and Acinetobacter baumannii . Recently, attention has also been focused on RND systems of other pathogenic bacteria, such as Burkholderia species , and more specifically, the Burkholderia cepacia complex (Bcc) . Bcc bacteria have been extensively investigated concerning their genome organization and the presence of several copies of MDR systems, especially the RND superfamily [20, 21].
Although rampant gene duplication played a key role in shaping this gene family, the tempo and mode of this process remains unknown. The genome of the model strain B. cenocepacia J2315 harbors several genes encoding RND proteins of different families [18, 20,21,22]. Two RND copies (RND 2 and 4) have been the focus of our attention since they provide a model to investigate the evolution of paralogs that are functionally distinct but with a high degree of sequence similarity. Previously, we have shown that the RND 2 coding gene is present only in some Bcc species, and a phylogenetic analysis revealed that the amino acid sequences of RND 2 and RND 4 proteins cluster together . Despite this high degree of sequence similarity, only the RND 4 protein plays a key role in the antibiotic resistance of B. cenocepacia J2315, since a deletion of the entire RND 4 operon reduces the intrinsic MDR antibiotic resistance of the mutant strain [23, 24]. Further, RND 4 appeared to be particularly important for antibiotic resistance of planktonic bacteria, but less relevant for resistance in biofilm bacteria . RND 4 also contributes to resistance to the biocide chlorhexidine  and the 2-thiopyridine anti-tubercular derivative, which is also effective against B. cenocepacia . This pump is also associated with the modulation of some virulence factors (motility, biofilm formation, chemotaxis and quorum sensing) . A proteomic analysis comparing the parental strain and the deletion mutant suggests a more general role for RND 4 in the physiology of B. cenocepacia cells . Conversely, the Burkholderia RND 2 protein is not expressed during growth on LB medium, and its role in antibiotic resistance can only be revealed when overexpressed in E. coli  .
In this study, we provide computational and experimental evidence supporting an evolutionary model on the functional diversification of these two RND superfamily members in Burkholderia spp., providing a framework to understand how these events could modulate antibiotic resistance.
Operon structure, comparative genomics and phylogeny
The RND 2 operon is located on B. cenocepacia J2315 chromosome 3. It consists of three genes organized in the following order: BCAS0766 (QU43_RS72495, 1239 bp), BCAS0765 (QU43_RS72490, 3192 bp), and BCAS0764 (QU43_RS72485, 1503 bp), encoding the MFP, RND permease, and OMP proteins, respectively. A LysR family transcriptional regulator (BCAS0767; QU43_RS72500, 881 bp) is located upstream of BCAS0766 and oriented in the same direction (Fig. 1a). Another gene, encoding an AraC family transcriptional regulator (BCAS0768; QU43_RS72510, 983 bp) is present upstream of the LysR regulator coding gene, but oriented in the opposite direction (Fig. 1a). The RND 4 operon is located on chromosome 1 and spans three genes: BCAL2822 (QU43_RS50725, 1275 bp), BCAL2821 (QU43_RS50720, 3201 bp), and BCAL2820 (QU43_RS50715, 1524 bp) encoding the MFP, RND permease, and OMP proteins, respectively. Thus, the RND 2 and RND 4 operons share identical gene organization. A gene encoding a TetR family transcriptional regulator (BCAL2823; QU43_RS50730, 641 bp), found upstream of BCAL2822, is oriented in the opposite direction (Fig. 1a). The degree of identity/similarity between the two operons is 92%, 87%, 95% of identity at the nucleotide level, with 90%, 85%, 96% identity and 93%, 91%, 96% similarity at the amino acid level when comparing the MFP, RND and OMP genes and proteins, respectively. The degree of DNA and amino acid sequence identity/similarity between the genes/proteins of the two operons (Fig. 1a) and the identical gene order strongly suggests these operons arose from a recent operon duplication event. No major evidence of purifying selection was found, as for almost all the pairwise comparisons among the RND, MFP, and OMP sequences, the dN/dS ratio was lower than 1 (Additional file 1 and Additional file 2).
A comparative genome analysis was performed on 797 Burkholderia genomes to evaluate the presence, distribution and conservation of RND 2 and RND 4 operons. RND 2 and 4 coding genes were found in 47 and 689 of the genomes analysed, respectively. Then the presence of the other two genes (mfp and omp) forming the operons was evaluated: the intact RND 2 operon occurred in 33 genomes, while the intact RND 4 operon was present in 472 (Additional file 3). The lack of the RND genes (or of the entire operons) in some strains could be due to the poor quality of the corresponding assemblies (draft genomes). In these genomes, operon fragments were found in contig boundaries. Therefore, incomplete operons were excluded from the downstream analysis. All genomes possessing the RND 2 operon were from Bcc bacterial species, but not all Bcc species harbor this operon. However, except for a few cases, the RND 2 operon co-occurred with RND 4 (Additional file 3). The phylogenetic tree of Fig. 1b (redrawn from the tree reported in De Smet et al. ) indicates the species in which a copy of the RND 2 operon was found. According to the phylogenetic distribution of RND4 and RND2 operons, the putative duplication event involving RND 4 and leading to RND 2 can be mapped inside the Bcc group (Fig. 1b).
The conservation of the genes encoding regulatory proteins associated with the RND 2 and RND 4 operons in B. cenocepacia J2315 was also investigated. A gene encoding a TetR family transcriptional regulator was always associated with the RND 4 operons identified (Additional file 3). LysR- and AraC-like transcriptional regulator genes present upstream of the J2315 RND 2 operon were also associated with all the RND 2 operons identified, and in the same order and orientation, except for B. diffusa (Fig. 1b) (Additional file 3). The RND 4 promoter regions were highly conserved among different strains and divergent from those of RND 2. The RND 2 promoter regions were also conserved among themselves, with the exception of B. diffusa (data not shown). In both cases, conserved motifs putatively recognized by regulatory proteins, were identified (Additional file 4).
Extra copies of the RND 2 operon can functionally replace the RND 4 operon
To assess the role of the two RND operons in antibiotic resistance and whether the RND 2 operon can functionally substitute RND 4, the B. cenocepacia J2315 RND 2 and RND 4 operons and their own promoter regions were cloned and inserted into pSCrhaB2 , giving rise to two recombinant plasmids (pEP_RND2_operon and pEP_RND4_operon, respectively) (Table 1). These plasmids were introduced in the parental J2315 strain and in a deletion mutant lacking the RND 4 operon (strain D4, Table 1). The antibiotic resistance profile of these strains was investigated by determining the Minimal Inhibitory Concentration (MIC) values of 14 antibiotics belonging to different classes. As a control, the MICs for a RND 2 operon deletion mutant were determined, and the same results were obtained as for J2315 (data not shown). The D4 strain was more sensitive than the parental J2315 to fluoroquinolones, chloramphenicol, tetracycline, rifampicin and novobiocin (Table 2), as previously reported [23, 24]. The other antibiotics tested are not substrates of RND 4 efflux pumps. As expected, the MIC values of D4(pEP_RND4_operon), in which the RND 4 pumps is restored, were identical to those of J2315 (except for levofloxacin, where the complementation of the RND 4 deletion was only partial). These results confirm that the reduction in MIC values was due to the absence of the RND 4 operon. The D4(pEP_RND2_operon) strain, which expresses a plasmid encoded RND 2 pump, showed an identical MIC profile as that of J2315 for ciprofloxacin and norfloxacin, a partial complementation for levofloxacin, sparfloxacin, tetracycline, rifampicin and novobiocin, and a MIC slightly higher for nalidixic acid and chloramphenicol. These results suggest that the RND 2 operon can complement the RND 4 deletion (although to a lesser extent for some antibiotics) despite not being able to carry out this function in the parental strain. The MIC values of J2315(pEP_RND2_operon) strain are not higher than those of J2315, while those of J2315(pEP_RND4_operon) are in some cases slightly higher, showing that additional copies of the two operons do not significantly increase the MIC values of the parental strain.
Native RND 2 is very poorly expressed in the D4 strain
The previous results indicate that extra copies of the RND 2 operon can functionally replace the RND 4 operon in the D4 deletion mutant by restoring resistance to some antibiotics. This raised the question of why the native RND 2 operon failed to complement the D4 mutation in vivo. To address this question, we investigated the expression of the RND 2 and RND 4 efflux protein genes by qRT-PCR in the presence of nalidixic acid, one of the common substrates of both pumps. In B. cenocepacia J2315, the estimated quantity of mRNA of RND 2 was about 100 times lower than that of RND 4 [3·10−7 and 3·10−5 normalized copies respectively (quantity mean normalized on the quantity mean of the 16S rRNA coding gene)], while no RND 4 operon expression was detected in the D4 mutant (Fig. 2). Further, the copy number of RND 2 mRNA remained low also in the D4 mutant strain (1·10−7 normalized copies), suggesting that the loss of RND 4 was not a condition sufficient to increase RND 2 expression. In addition, a strong increase of RND 2 expression was detected in the strain D4(pEP_RND2_operon) (2·10−5 normalized copies) (Fig. 2), consistent with the finding that this strain re-acquired resistance to some of the antibiotics tested. The extra-copies of the RND 2 operon in J2315 increased the expression of both RND 4 (6·10−4 normalized copies) and RND 2 (4·10−5 normalized copies); apparently, this should result in a higher degree of resistance to antibiotics. However, the strong increase of RND 2 expression was masked by the presence of RND 4 efflux pump, which was the main one responsible for the resistance to the antibiotics tested. The addition of extra-copies of RND 4 to J2315 did not significantly increase the expression of RND 4 (1·10−4 normalized copies) nor RND 2 (6·10−7 normalized copies).
Spontaneous mutants with high antibiotic resistance can be obtained from the D4 strain through directed evolution experiments
The ability of RND 2 to partially replace the activity of the RND 4 operon through an increase in expression suggests that under selective pressure it should be possible to isolate spontaneous mutants of strain D4 with MIC values similar to those of strain D4(pEP_RND2_operon). Therefore, we set up a directed evolution experiment with the B. cenocepacia D4 strain to obtain spontaneous mutants resistant to sublethal concentrations of chloramphenicol (one of the antibiotics that is a substrate of both RND pumps). Two spontaneous mutants, D4/C18 and D4/C20, were obtained after 15 days in the presence of chloramphenicol. These mutants were examined for antimicrobial resistance and the MIC values against different antibiotics for the two mutants were similar (Table 2). Increased resistance against chloramphenicol and novobiocin compared to both D4 and J2315 was observed. For quinolones/fluoroquinolones the behavior was different depending on the antibiotic: nalidixic acid, ciprofloxacin, levofloxacin and norfloxacin had MIC values higher than that of original D4 strain and similar to those of J2315 strain, while sparfloxacin resistance remained that of D4 strain. Therefore, D4/C18 and D4/C20 showed increased multiple antibiotic resistance relative to D4 and not only resistance to the antibiotic used for selection, suggesting the activation of an efflux mechanism.
An insertion activates RND 2 expression in the D4/C18 and D4/C20 strains
RAPD fingerprinting analysis confirmed that the D4/C18 and D4/C20 mutants had a profile identical to D4, confirming their clonality (data not shown). We then extracted and sequenced the genomic DNA of strains D4, D4/C18, and D4/C20 to identify the mutation(s) arising under selective pressure (Additional file 5). The genome sequences of the two mutants were compared to those of D4 and J2315 for the presence of SNPs and/or large insertions/deletions. No SNPs were found in the two mutants at the whole genome level, but several deletions and insertions of different size were detected (Table 3). The only mutation common to the D4/C18 and D4/C20 strains, associated with the common phenotype, is a 1300-bp insertion at position 867,051 of chromosome 3, confirmed by PCR and DNA sequencing. This fragment carries predicted transposase and integrase genes both present on chromosome 1 (BCAL2755 = QU43_RS50390 and BCAL2756 = QU43_RS50395) and chromosome 2 (BCAM1929 = QU43_RS63980 and BCAM1930 = QU43_RS63985) of J2315. The insertion is located 116 bp upstream of the araC family transcriptional regulator gene (BCAS0768 = QU43_RS72510) and oriented in the same direction, and 482 bp upstream of the lysR family transcriptional regulator gene (BCAS0767 = QU43_RS72500), oriented divergently. These two genes encode regulatory proteins that are associated in the genome with the RND 2 operon. Such genes could be involved in its regulation of expression. To check the possible effects of this insertion on the surrounding genes, the expression levels of both genes and of the RND 2 coding gene were evaluated in D4, D4/C18 and D4/C20 strains (Fig. 3). In both D4/C18 and D4/C20 the expression of the araC regulator gene was double that of D4, while the expression of both the lysR regulator and RND 2 genes increased significantly (t-test, p < 0.05). We concluded that the increased resistance observed in D4/C18 and D4/C20 was due to increased RND 2 operon expression.
We investigated how paralog RND operons could have evolved and functionally diversified, and how these events could influence the onset of novel resistance phenotypes. The search for RND 2 and RND 4 operons in 797 Burkholderia genomes revealed that the RND 4 operon is highly conserved, being present in almost all the genomes examined. The high conservation of RND 4 is consistent with the critical role of this operon in multiple antibiotic resistance, virulence, and other cellular processes [23,24,25,26, 28]. Further, a gene encoding a TetR family transcriptional regulator located upstream of RND 4 is highly conserved, suggesting that regulatory elements controlling the RND 4 operon expression are also common to all Burkholderia species.
In contrast, the RND 2 operon is only found in a few Bcc species and likely resulted from a duplication event of the entire RND 4 operon. The conclusion that RND 2 originated by a duplication of RND 4 and not by duplication of other RND operons sharing the very same gene organization  is supported by the high degree of sequence identity between RND 2 and RND 4, which is much higher than those existing between RND 2 or RND 4 and the other RND operons. Lastly, the phylogenetic distribution of RND 2 operon allowed us to map the putative duplication event inside the Bcc group (Fig. 1b).
The RND 2 operon has retained most of the ancestral substrate specificity, in that it maintained the ability to extrude the same antimicrobial compounds extruded by RND 4. This hypothesis is supported also by the ability of RND 2 to restore the resistance to some antibiotics in the D4 mutant. However, the recovery of resistance required overexpression of RND 2, that can be obtained either by introducing extra-copies of RND 2 or by mutations altering its regulation. Confirming this scenario, qRT-PCR experiments revealed that RND 2 is expressed at very low level in LB medium even in the presence of a possible inducer, in both B. cenocepacia J2315 and the D4 mutant. The need for extra-copies of the RND 2 operon in the D4 mutant or of a mutation to increase the expression of the RND 2 operon suggests that antibiotics are not the stimulus that leads to the activation RND 2. According to these hypothesis, different expression “circuits” should control the two operons. From an evolutionary viewpoint, this means that the duplication of the RND 4 operon also involved parallel or subsequent genetic rearrangements resulting in new regulatory elements (e.g. promoters) located upstream of the MFP gene of the RND 2 operon (subfunctionalization) (Fig. 4). This new regulatory control may have involved the appearance of two new genes, encoding AraC-like and LysR-type regulators, respectively (Fig. 1a). These two genes are present in all genomes harboring an RND 2 operon, except for B. diffusa (Fig. 1b). This suggests that the molecular rearrangements that led to the localization of these two genes close to the RND2 operon, might have occurred after the duplication event, since B. diffusa is phylogenetically located at the base of the cluster in which the duplication occurred (Fig. 1b). The origin of these two genes is still unclear. Searches carried out in all the available prokaryotic genomes did not retrieve any sequence with a degree of identity/similarity sufficiently high to suggest a common origin. Regardless, this regulatory circuit, once acquired, has been maintained, since the two regulatory genes are present and in the same order in all genomes harboring an entire RND 2 operon.
Expression data shed light on the regulation of the RND 2. In the presence of extra copies of RND 2, supplied by the recombinant plasmid pEP_RND2_operon, the expression of RND 2 in J2315 (pEP_RND2_operon) strain is higher than that in the parental strain. A plausible scenario could be that a constitutively expressed repressor controls the RND 2 operon. We can then hypothesize that when extra copies of the RND 2 operon are present, the repressor concentration might be not sufficient to block RND 2 expression. Accordingly, mutations increasing the expression of a regulator gene would have similar effects. This hypothesis was supported by the results from the directed evolution experiments. In fact, increased RND 2 expression was parallel to a strong increase of the associated lysR gene expression. Also, the gene encoding the AraC-type protein might play a role in the RND 2 expression, since the expression of this gene increases in both mutants D4/C18 and D4/C20.
In summary, our results show that the Bcc RND 2 operon most likely arose from a duplication of the RND 4 operon, but although RND 2 retains almost the same functions in antibiotics extrusion as RND 4, its expression in Burkholderia cells is much lower. Therefore, a different local mechanism regulating the expression of the RND 2 operon suggests that changes in regulatory circuits (rather than functional mutations) may influence the role of different RND paralogs. Interestingly, the fact that the RND 2 operon is encoded by B. cenocepacia J2315 chromosome 3 further supports the subfunctionalization hypothesis concerning the evolution of these two efflux systems. Chromosome 3 is a megaplasmid encoding several critical features for pathogenicity and chronic infection, but dispensable for bacterial viability in vitro [31, 32]. This agrees with the notion that the RND 2 operon is not completely diversified from its homolog operon (i.e. RND 4), and consequently, it is not yet essential for B. cenocepacia cells. It is plausible that the RND 2 operon may be required only under certain conditions related to chronic infection. Our directed evolution experiment demonstrated that under selective pressure, mutations may arise that rewire the local regulatory circuit of the RND 2 operon, allowing it to respond to the presence of antibiotics. A similar mechanism was identified in E. coli cells responding to osmotic stress , suggesting that expression rewiring is an evolutionary answer to stress conditions in bacteria.
Our study contributes to an explanation of how bacterial regulatory networks can evolve to allow partial functional diversification of paralogous genes, thus becoming key players in maintaining a high genomic (and consequently phenotypic) plasticity. Understanding the regulatory network underpinning the expression of the RND multimeric complexes in Bcc can add new information explaining the evolution of drug resistant phenotypes and also potentially new strategies for their control.
Bacterial strains and growth conditions
The bacterial strains and plasmids used in this work are listed in Table 1. Bacteria were grown under aerobic condition at 37 °C in Luria-Bertani (LB) agar or broth. Antibiotic concentration used were 100 μg/ml ampicillin, 50 μg/ml and 50 or 100 μg/ml trimethoprim for E. coli and B. cenocepacia, respectively. Both antibiotics were purchased from Sigma-Aldrich S.r.l.
The genome sequence of B. cenocepacia J2315 available from GenBank database (GCA_000009485.1) , was used for primers design. Due to their high size and GC content, the two RND operons, RND 2 (old locus tags: BCAS0764–BCAS0766, new locus tags: QU43_RS72485, QU43_RS72490, QU43_RS72495) and RND 4 (old locus tags: BCAL2820–BCAL2822, new locus tags: QU43_RS50715, QU43_RS50720, QU43_RS50725), were cloned using a two-step strategy. Firstly, a unique restriction site (BamHI for operon RND 2 and KpnI for operon RND 4) was identified in the sequences of the two operons. In this way, the two operons were split into two parts, and two divergent primers (R2_1/F2_2 and R4_1/F4_2 for operon RND 2 and RND 4 respectively, Additional file 6) were designed in correspondence of these two restriction sites. Then, for each operon a forward primer for the amplification of the first part of the operon (including the putative promoter region) and a reverse primer for the amplification of the second part of the operon (including the putative transcription terminator regions) were designed (Additional file 6). These primers contained the restriction sites necessary for the subsequent recovery of the operon (in both cases an NdeI and HindIII restriction site were introduced in the forward and reverse primers, respectively). The primer pair F2_1/R2_1 was used to amplify a fragment of 3479 bp containing the putative promoter region of operon RND 2, the gene coding for the MFP (BCAS0766, QU43_RS72495) and a portion of the RND protein-encoding gene (BCAS0765, QU43_RS72490). The primer pair F2_2/R2_2 allowed for the amplification of a 2699 bp fragment containing the second half of the RND coding gene (BCAS0765, QU43_RS2490) and the OMP coding gene (BCAS0764, QU43_RS72485). For the RND 4 operon, the primer pair F4_1/R4_1 was used to amplify a 1937 bp fragment containing the putative promoter region, the MFP coding gene (BCAL2822, QU43_RS50725) and the first half of the RND protein-encoding gene (BCAL2821, QU43_RS50720). The primer pair F4_2/R4_2 for the amplification to a 4351 bp fragment containing the second half of the RND coding gene (BCAL2821, QU43_RS50720) and the OMP coding gene (BCAL2820, QU43_RS50715).
PCR amplification of each region was performed in a 50 μL reaction mixture containing an aliquot of 2 μL of purified DNA [prepared employing the NucleoSpin® Tissue extraction kit (MACHEREY-NAGEL GmbH & Co. KG) following the manufacturer’s instructions], 10 μL of 5X reaction buffer (5X Phusion HF Buffer, Thermo Fisher Scientific), 0.5% DMSO, 0.5 μM of each primer, 200 μM dNTPs and 1 U of Taq DNA polymerase (Phusion High-Fidelity DNA Polymerase, Thermo Fisher Scientific). In order to reduce the number of non-specific PCR products, an amplification program with several increasing annealing temperatures was used . For the two primer pairs F2_1/R2_1 and F2_2/R2_2, cycle condition were 98 °C for 3′, followed by 30 cycles of 98 °C for 45″, 60 °C for 45″, 62.5 °C for 45″, 65 °C for 45″, 67.5 °C for 45″, 70 °C for 45″, 72 °C for 2′, and a final extension of 10′ at 72 °C. For the primer pair F4_1/R4_1 cycle condition were 98 °C for 3′, followed by 30 cycles of 98 °C for 45″, 67.5 °C for 45″, 70 °C for 45″, 72 °C for 1′, and a final extension of 10′ at 72 °C. For the primer pair F4_2/R4_2 cycling conditions were 98 °C for 3′, followed by 30 cycles of 98 °C for 45″, 65 °C for 45″, 72 °C for 2′ and 30″, and a final extension of 10′ at 72 °C. Amplification products were loaded on 0.6% agarose gels and amplicons were excised from agarose gel and purified using the MinElute gel extraction kit (Qiagen) according to the manufacturer’s instructions.
Then, 3′ A-overhangs were added to the four purified PCR products, in order to permit their insertion into the appropriate plasmid cloning vector (see below). To this purpose, 0.2 mM dATP, 1X reaction buffer (10X DreamTaq Buffer, Thermo Fisher Scientific) and 2 U of Taq polymerase (DreamTaq DNA Polymerase, Thermo Fisher Scientific) were added to the purified DNA in a final volume of 10 μl, and incubated for 20′ at 72 °C. After that, the four amplicons were cloned into the pGEM®-T Easy Vector (Promega), using the manufacturer’s instructions. The obtained plasmids (Additional file 7) were introduced into E. coli DH5α by electroporation . The integrity of the cloned fragments was verified by sequencing of the plasmids at the Genechron laboratory (Ylichron Srl, Italy) using the primers reported in Additional file 6.
The entire RND 2 and RND 4 operons were then assembled in the pSCrhaB2 vector . First the two plasmids pGEM_RND2_second_part and pGEM_RND4_second_part (Additional file 7) were digested with BamHI-HindIII and KpnI-HindIII, respectively. The obtained fragments, corresponding to the second parts of the two operons, were ligated into the BamHI-HindIII and KpnI-HindIII digested pSCrhaB2, respectively, to yield the pEP_RND2_second_part and pEP_RND4_second_part plasmids (Additional file 7). Then, the two plasmids pGEM_RND2_first_part and pGEM_RND4_first_part (Additional file 7) were digested with NdeI-BamHI and NdeI-KpnI, respectively, and the obtained fragments were ligated into the NdeI-BamHI and NdeI-KpnI digested pEP_RND2_second_part and pEP_RND4_second_part, respectively, to obtain the pEP_RND2_operon and pEP_RND4_operon plasmids (Table 1 and Additional file 7). All the restriction enzymes were purchased from Fermentas (Thermo Fisher Scientific), while for the ligation reactions the T4-DNA ligase included in the pGEM®-T Easy Vector kit (Promega) was used.
The two plasmids pEP_RND2_operon and pEP_RND4_operon were introduced in E. coli SY327  by electroporation  and then mobilized into B. cenocepacia J2315 by triparental mating, using the helper plasmid pRK2013 [37, 38]. The presence of plasmids in the final strains was controlled.
The MIC of different classes of antibiotics were determined: fluoroquinolones (ciprofloxacin, levofloxacin, norfloxacin, sparfloxacin, nalidixic acid), aminoglycosides (kanamycin), ansamycins (streptomycin), β-lactams (ampicillin), macrolides (erythromycin), tetracycline, rifampicin, chloramphenicol, gramidicin and novobiocin. All the antibiotics were purchased from Sigma-Aldrich S.r.l..
The MIC determination protocol was adapted from Ulrich et al. . Briefly, an Over-Night (ON) culture in LB broth (and trimethoprim 50 μg/ml for strains carrying plasmids) of each strain was diluted to a concentration of about 106 CFU/ml. In each well of a 96-well plate, 50 μl of these bacterial suspensions were added to plate wells containing LB medium with the appropriate concentration of antibiotics (and trimethoprim at a final concentration of 50 μg/ml for strains carrying plasmids) to obtain a final bacterial concentration of about 5 × 104 CFU/well. Plates were incubated for 48 h at 37 °C statically and MICs determined both visually and based on the OD600 using a Tecan Infinite M200 plate reader (Tecan, San Jose, CA). All MICs were determined in triplicate, and the MIC was defined as the lowest concentration of antibiotics that prevented any detectable growth.
RNA extraction and reverse transcription
B. cenocepacia strains J2315, D4, D4(pEP_RND2_operon), D4(pEP_RND4_operon), J2315(pEP_RND2_operon), and J2315(pEP_RND4_operon) were grown in LB broth in the presence of 4 μg/ml of nalidixic acid (half the MIC of the D4 strain), and trimethoprim 50 μg/ml for strains carrying plasmids, until an OD600 ~ 0.5–0.6. 500 μl of each culture were treated with the RNA protect bacteria reagent (Qiagen) and total RNA was extracted using the RNeasy Mini Kit (Qiagen), using the manufacturer’s instruction. DNA was removed from the sample using the RNase-free DNase set (Qiagen).
Extracted RNA was reverse-transcribed using the Superscript II Reverse Transcriptase (Invitrogen) and Random primers (Invitrogen) following the manufacturer’s instruction.
B. cenocepacia D4, D4/C18 and D4/C20 strains were grown in LB broth until an OD600 ~ 0.5–0.6 and RNA was extracted and reverse-transcribed from 500 μl of each of these cultures following the same procedure.
Quantitative real-time PCR (qRT-PCR)
qRT-PCR reactions were performed in 10 μl reaction mixture containing an aliquot of 1 μl of cDNA, 5 μl of Maxima Syber Green/Rox qPCR Master Mix (2X) (ThermoScientific) and 1 μg/ml of each primer. For B. cenocepacia J2315, D4, D4(pEP_RND2_operon), D4(pEP_RND4_operon), J2315(pEP_RND2_operon), and J2315(pEP_RND4_operon), each cDNA, amplification was carried out using both the primer pairs RND_2_for/RND_2_rev and RND_4_for/RND_ 4_rev (Additional file 6), while 16S RNA encoding gene was used as reference. For the B. cenocepacia D4, D4/C18 and D4/C20 strains amplifications were performed using the three primer pairs RND_2_for/RND_2_rev, BCAS0767_for/BCAS0767_rev and BCAS0768_for/BCAS0768_rev (Additional file 6). The expression of each gene in the mutant strain was compared with the expression of the same gene in the D4 strain; 16S RNA encoding gene was used as reference.
Each sample was spotted in triplicate and known amounts of DNA of the B. cenocepacia J2315 strain (1–0.1-0.01-0.001 ng) were added to obtain a standard curve for each primer pairs. The reactions were performed on a QuantStudio™ 7 Flex Real-Time PCR System (Applied Biosystems by Life Technologies). Cycling conditions were: hold stage [50 °C for 2′ and 95 °C for 10′], PCR stage [40 cycles of: 95 °C for 30″, 57 °C or 60 °C for 30″ (57 °C for the primer pairs RND_2_for/RND_2_rev and RND_4_for/RND_ 4_rev, 60 °C for other pairs), 72 °C for 15″], melt curve stage [95 °C for 15″ and 60 °C for 1′].
Directed evolution experiment and mutants fingerprinting
B. cenocepacia D4 strain was grown ON (16 h) in LB broth with shaking. 100 μl of this culture (~ 108–109 cells) and of a 10−1 dilution were plated on LB agar containing different concentration of chloramphenicol (16, 32, 64, 128 μg/ml). All plates were checked every day, for 15 days, and all of the colonies that grew, were further selected on plates with 128 μg/ml chloramphenicol. In this way 2 spontaneous resistant mutants selected on chloramphenicol (named D4/C18 and D4/C20, Table 1) were obtained.
DNA extraction and genome sequencing
Genomic DNA of B. cenocepacia D4, D4/C18 and D4/C20 was extracted using the CTAB protocol previously described in Perrin et al. . Whole genome shot-gun sequencing was performed by the Institute of Applied Genomics and IGA Technology Services S.r.l. (University of Udine, Italy) using an Illumina (Solexa) HiSeq2500.
The genome sequences determined are available in GenBank, accession numbers: B. cenocepacia D4 (SRR3736982), B. cenocepacia D4/C18 (SRR3737008), B. cenocepacia D4/C20 (SRR3737019).
Deletions and insertions were confirmed by PCR amplification of the surrounding regions of the hypothetical mutation using the primers reported in Additional file 6. PCR amplifications were performed using the DreamTaq DNA Polymerase (Thermo Fisher Scientific) in a 10 μl reaction mixture containing 1 μl of same purified DNA, 1X reaction buffer, 0.5 μM of each primer, 200 μM dNTPs and 2 U of Taq polymerase. Cycling condition were: 95 °C for 3′, followed by 30 cycles of 95 °C for 10″, 60 °C for 30″, 72 °C for xx” (depending on the size), and a final extension of 10′ at 72 °C.
Amplification products were purified as reported above and sequenced at the Genechron laboratory (Ylichron S.r.l., Italy).
Blast search of RND 2 and 4 operons
All the genomic sequences of Burkholderia representatives available from GenBank (up to October 25th, 2016) were downloaded using in-house bash scripts: 797 sequences were obtained. When not present in GenBank, genome annotation was obtained by applying the PROKKA  annotation software to the genome sequence.
The amino acid sequences of RND 4 (gi: 206,561,158) and RND 2 (gi: 197,295,595) from B. cenocepacia J2315 were used as queries for a blastp search in all the accessions obtained. For each paralog, only best hits with e-values <= 1e−15 were considered. In this way, RND 2 and 4 coding genes were found in 47 and 689 of the genomes analyzed, respectively. Neighboring coding DNA sequences (CDSs) were then compared with the amino acid sequences of OMP and MFP from B. cenocepacia J2315 (using blastp, considering only the best hits with e-value <= 1e−15) to identify the corresponding homologs in the target genome. The entire RND 2 operon was found in 33 genomes, while the entire RND 4 operon in 472 (Additional file 3). The homologs of AraC (gi: 197,295,598), TetR (gi: 206,561,160) and LysR (gi: 197,295,597) were found in a similar way.
Sequencing reads analysis and variants calling
Reads quality was assessed using FastQC toolkit (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Reads trimming was performed using the dynamic trimming approach implemented in SolexQA , using a Phred score of 20 as the base call quality threshold. Trimmed reads were then mapped on the reference genome using BWA with default parameters. Only those SNPs that (i) were not present in the WT strain, (ii) were supported by, at least, five reads, and (iii) had a support greater than 50% were considered for further analysis. SNPs calling was performed using VarScan 2.4.0 . Large variants (insertion and deletion) calling was performed on the sorted BAM file produced by BWA using Pilon 1.16 . All SAM and BAM files manipulation was performed using SAM tools .
Burkholderia cepacia complex
- MDR EP:
Multidrug resistance efflux pump
Multidrug antibiotic resistance
Membrane fusion protein
Minimal Inhibitory Concentration
Next generation sequencing
Outer membrane channel protein
Quantitative Real Time PCR
W.H.O. Antimicrobial resistance: global report on surveillance. Geneva, Switzerland: World Health Organization; 2014.
Venter H, Mowla R, Ohene-Agyei T, Ma S. RND-type drug e ffl ux pumps from Gram-negative bacteria: molecular mechanism and inhibition. Front Microbiol. 2015;6:377.
Li XZ, Plesiat P, Nikaido H. The challenge of efflux-mediated antibiotic resistance in Gram-negative bacteria. Clin Microbiol Rev. 2015;28(2):337–418.
Muller B, Borrell S, Rose G, Gagneux S. The heterogeneous evolution of multidrug-resistant mycobacterium tuberculosis. Trends in genetics : TIG. 2013;29(3):160–9.
Martinez JL, Sanchez MB, Martinez-Solano L, Hernandez A, Garmendia L, Fajardo A, Alvarez-Ortega C. Functional role of bacterial multidrug efflux pumps in microbial natural ecosystems. FEMS Microbiol Rev. 2009;33(2):430–49.
Serres MH, Kerr AR, McCormack TJ, Riley M. Evolution by leaps: gene duplication in bacteria. Biol Direct. 2009;4:46.
Saier MH Jr, Paulsen IT, Sliwinski MK, Pao SS, Skurray RA, Nikaido H. Evolutionary origins of multidrug and drug-specific efflux pumps in bacteria. FASEB J. 1998;12(3):265–74.
Kondrashov FA, Kondrashov AS. Role of selection in fixation of gene duplications. J Theor Biol. 2006;239(2):141–51.
Craven SH, Neidle EL. Double trouble: medical implications of genetic duplication and amplification in bacteria. Future Microbiol. 2007;2(3):309–21.
Sandegren L, Andersson DI. Bacterial gene amplification: implications for the evolution of antibiotic resistance. Nat Rev Microbiol. 2009;7(8):578–88.
Kondrashov FA. Gene duplication as a mechanism of genomic adaptation to a changing environment. Proc Biol Sci. 2012;279(1749):5048–57.
Lynch M, Force A. The probability of duplicate gene preservation by subfunctionalization. Genetics. 2000;154(1):459–73.
Force A, Cresko WA, Pickett FB, Proulx SR, Amemiya C, Lynch M. The origin of subfunctions and modular gene regulation. Genetics. 2005;170(1):433–46.
Alvarez-Ortega C, Olivares J, Martinez JL. RND multidrug efflux pumps: what are they good for? Front Microbiol. 2013;4:7.
Martinez JL. Antibiotics and antibiotic resistance genes in natural environments. Science. 2008;321(5887):365–7.
Saier MH Jr, Paulsen IT. Phylogeny of multidrug transporters. Semin Cell Dev Biol. 2001;12(3):205–13.
Poole K, Krebes K, McNally C, Neshat S. Multiple antibiotic resistance in Pseudomonas aeruginosa: evidence for involvement of an efflux operon. J Bacteriol. 1993;175(22):7363–72.
Podnecky NL, Rhodes KA, Schweizer HP. E ffl ux pump-mediated drug resistance in Burkholderia. Front Microbiol. 2015;6:305.
De Smet B, Mayo M, Peeters C, Zlosnik JE, Spilker T, Hird TJ, LiPuma JJ, Kidd TJ, Kaestli M, Ginther JL et al: Burkholderia stagnalis sp. nov. and Burkholderia territorii sp. nov., two novel Burkholderia cepacia complex species from environmental and human sources. Int J Syst Evol Microbiol 2015, 65(7):2265-2271.
Perrin E, Fondi M, Papaleo MC, Maida I, Buroni S, Pasca MR, Riccardi G, Fani R. Exploring the HME and HAE1 efflux systems in the genus Burkholderia. BMC Evol Biol. 2010;10:164.
Perrin E, Fondi M, Papaleo MC, Maida I, Emiliani G, Buroni S, Pasca MR, Riccardi G, Fani R. A census of RND superfamily proteins in the Burkholderia genus. Future Microbiol. 2013;8(7):923–37.
Holden MT, Seth-Smith HM, Crossman LC, Sebaihia M, Bentley SD, Cerdeno-Tarraga AM, Thomson NR, Bason N, Quail MA, Sharp S, et al. The genome of Burkholderia cenocepacia J2315, an epidemic pathogen of cystic fibrosis patients. J Bacteriol. 2009;191(1):261–77.
Buroni S, Pasca MR, Flannagan RS, Bazzini S, Milano A, Bertani I, Venturi V, Valvano MA, Riccardi G. Assessment of three resistance-nodulation-cell division drug efflux transporters of Burkholderia cenocepacia in intrinsic antibiotic resistance. BMC Microbiol. 2009;9:200.
Bazzini S, Udine C, Sass A, Pasca MR, Longo F, Emiliani G, Fondi M, Perrin E, Decorosi F, Viti C, et al. Deciphering the role of RND efflux transporters in Burkholderia cenocepacia. PLoS One. 2011;6(4):e18902.
Buroni S, Matthijs N, Spadaro F, Van Acker H, Scoffone VC, Pasca MR, Riccardi G, Coenye T: Differential roles of RND efflux pumps in antimicrobial drug resistance of sessile and planktonic Burkholderia cenocepacia cells. Antimicrob Agents Chemother 2014, 58(12):7424-7429.
Coenye T, Van Acker H, Peeters E, Sass A, Buroni S, Riccardi G, Mahenthiralingam E: Molecular mechanisms of chlorhexidine tolerance in Burkholderia cenocepacia biofilms. Antimicrob Agents Chemother 2011, 55(5):1912-1919.
Scoffone VC, Spadaro F, Udine C, Makarov V, Fondi M, Fani R, De Rossi E, Riccardi G, Buroni S: Mechanism of resistance to an antitubercular 2-thiopyridine derivative that is also active against Burkholderia cenocepacia. Antimicrob Agents Chemother 2014, 58(4):2415-2417.
Gamberi T, Rocchiccioli S, Papaleo MC, Magherini F, Citti L, Buroni S, Bazzini S, Udine C, Perrin E, Modesti A, et al. RND-4 efflux transporter gene deletion in Burkholderia cenocepacia J2315: a proteomic analysis. J Proteome Sci Comput Biol 2013. 2013;2:1.
Guglierame P, Pasca MR, De Rossi E, Buroni S, Arrigo P, Manina G, Riccardi G: Efflux pump genes of the resistance-nodulation-division family in Burkholderia cenocepacia genome. BMC Microbiol 2006, 6:66.
Cardona ST, Valvano MA. An expression vector containing a rhamnose-inducible promoter provides tightly regulated gene expression in Burkholderia cenocepacia. Plasmid. 2005;54(3):219–28.
Agnoli K, Schwager S, Uehlinger S, Vergunst A, Viteri DF, Nguyen DT, Sokol PA, Carlier A, Eberl L. Exposing the third chromosome of Burkholderia cepacia complex strains as a virulence plasmid. Mol Microbiol. 2012;83(2):362–78.
Agnoli K, Frauenknecht C, Freitag R, Schwager S, Jenul C, Vergunst A, Carlier A, Eberl L. The third replicon of members of the Burkholderia cepacia complex, plasmid pC3, plays a role in stress tolerance. Appl Environ Microbiol. 2014;80(4):1340–8.
Stoebel DM, Hokamp K, Last MS, Dorman CJ. Compensatory evolution of gene regulation in response to stress by Escherichia Coli lacking RpoS. PLoS Genet. 2009;5(10):e1000671.
Choi KH, Kumar A, Schweizer HP. A 10-min method for preparation of highly electrocompetent Pseudomonas Aeruginosa cells: application for DNA fragment transfer between chromosomes and plasmid transformation. J Microbiol Methods. 2006;64(3):391–7.
Sambrook J, Fritsch EF, Maniatis T. Molecular cloning : a laboratory manual. 2nd ed. Cold Spring Harbor, N.Y: Cold Spring Harbor Laboratory; 1989.
Miller VL, Mekalanos JJ. A novel suicide vector and its use in construction of insertion mutations: osmoregulation of outer membrane proteins and virulence determinants in Vibrio cholerae requires toxR. J Bacteriol. 1988;170(6):2575–83.
Figurski DH, Helinski DR. Replication of an origin-containing derivative of plasmid RK2 dependent on a plasmid function provided in trans. Proc Natl Acad Sci U S A. 1979;76(4):1648–52.
Craig FF, Coote JG, Parton R, Freer JH, Gilmour NJ. A plasmid which can be transferred between Escherichia Coli and Pasteurella haemolytica by electroporation and conjugation. J Gen Microbiol. 1989;135(11):2885–90.
Ulrich RL, Deshazer D, Kenny TA, Ulrich MP, Moravusova A, Opperman T, Bavari S, Bowlin TL, Moir DT, Panchal RG. Characterization of the Burkholderia thailandensis SOS response by using whole-transcriptome shotgun sequencing. Appl Environ Microbiol. 2013;79(19):5830–43.
Mocali S, Galeffi C, Perrin E, Florio A, Migliore M, Canganella F, Bianconi G, Di Mattia E, Dell'Abate MT, Fani R et al: Alteration of bacterial communities and organic matter in microbial fuel cells (MFCs) supplied with soil and organic fertilizer. Appl Microbiol Biotechnol 2013, 97(3):1299-1315.
Mori E, Lio P, Daly S, Damiani G, Perito B, Fani R. Molecular nature of RAPD markers from Haemophilus influenzae rd genome. Res Microbiol. 1999;150(2):83–93.
Di Cello F, Bevivino A, Chiarini L, Fani R, Paffetti D, Tabacchioni S, Dalmastri C: Biodiversity of a Burkholderia cepacia population isolated from the maize rhizosphere at different plant growth stages. Appl Environ Microbiol 1997, 63(11):4485-4493.
Perrin E, Fondi M, Maida I, Mengoni A, Chiellini C, Mocali S, Cocchi P, Campana S, Taccetti G, Vaneechoutte M, et al. Genomes analysis and bacteria identification: the use of overlapping genes as molecular markers. J Microbiol Methods. 2015;117:108–12.
Seemann T. Prokka: rapid prokaryotic genome annotation. Bioinformatics. 2014;30(14):2068–9.
Cox MP, Peterson DA, Biggs PJ. SolexaQA: at-a-glance quality assessment of Illumina second-generation sequencing data. BMC bioinformatics. 2010;11:485.
Koboldt DC, Zhang Q, Larson DE, Shen D, McLellan MD, Lin L, Miller CA, Mardis ER, Ding L, Wilson RK. VarScan 2: somatic mutation and copy number alteration discovery in cancer by exome sequencing. Genome Res. 2012;22(3):568–76.
Walker BJ, Abeel T, Shea T, Priest M, Abouelliel A, Sakthikumar S, Cuomo CA, Zeng Q, Wortman J, Young SK, et al. Pilon: an integrated tool for comprehensive microbial variant detection and genome assembly improvement. PLoS One. 2014;9(11):e112963.
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R, Genome project data processing S. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25(16):2078–9.
The authors are grateful to George C diCenzo for the mother tongue language revision of the manuscript.
This work was supported by the Ente Cassa di Risparmio di Firenze [project 2016.0936]. Elena Perrin was founded by a “Fondazione Adriano Buzzati-Traverso” fellowship.
Availability of data and materials
The genome sequence of B. cenocepacia J2315 used in this work is available in GenBank database (GCA_000009485.1). The genome sequences determined in this work are available in GenBank database with the following accession number: B. cenocepacia D4 (SRR3736982), B. cenocepacia D4/C18 (SRR3737008), B. cenocepacia D4/C20 (SRR3737019).
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Calculation of synonymous/non-synonymous substitution. (ZIP 34380 kb)
Materials and Methods of the calculation of synonymous/non-synonymous substitution and of the analysis of the promoter regions. (DOC 26 kb)
Distribution of the RND 2 and 4 operons and of their regulator in Burkholderia genomes. (XLS 92 kb)
Analysis of promoter regions. (ZIP 1251 kb)
Quality and depth of the sequencing run on B. cenocepacia mutants. (XLS 29 kb)
Primers used in this work. (XLS 40 kb)
Complete list of all the plasmids and strains used in this work. (XLS 35 kb)
About this article
Cite this article
Perrin, E., Fondi, M., Bosi, E. et al. Subfunctionalization influences the expansion of bacterial multidrug antibiotic resistance. BMC Genomics 18, 834 (2017). https://doi.org/10.1186/s12864-017-4222-4