Modulation of Global Transcriptional Regulatory Networks as a Strategy for Increasing Kanamycin Resistance of the Translational Elongation Factor-G Mutants in Escherichia coli

Evolve and resequence experiments have provided us a tool to understand bacterial adaptation to antibiotics. In our previous work, we used short-term evolution to isolate mutants resistant to the ribosome targeting antibiotic kanamycin, and reported that Escherichia coli develops low cost resistance to kanamycin via different point mutations in the translation Elongation Factor-G (EF-G). Furthermore, we had shown that the resistance of EF-G mutants could be increased by second site mutations in the genes rpoD/cpxA/topA/cyaA. Mutations in three of these genes had been discovered in earlier screens for aminoglycoside resistance. In this work, we expand our understanding of these second site mutations, the goal being to understand how these mutations affect the activities of the mutated gene products to confer resistance. We show that the mutation in cpxA most likely results in an active Cpx stress response. Further evolution of an EF-G mutant in a higher concentration of kanamycin than what was used in our previous experiments identified the cpxA locus as a primary target for a significant increase in resistance. The mutation in cyaA results in a loss of catalytic activity and probably results in resistance via altered CRP function. Despite a reduction in cAMP levels, the CyaAN600Y mutant has a transcriptome indicative of increased CRP activity, pointing to an unknown role for CyaA and / or cAMP in gene expression. From the transcriptomes of double and single mutants, we describe the epistasis between the mutation in EF-G and these second site mutations. We show that the large scale transcriptomic changes in the topoisomerase I (FusAA608E-TopAS180L) mutant likely result from increased negative supercoiling in the cell. Finally, genes with known roles in aminoglycoside resistance were present among the misregulated genes in the mutants.


Identification of SNPs from deep-sequencing data
The following pipeline was used to identify SNPs from genomic-DNA / RNA deep sequencing data. Reads were trimmed using CUTADAPT (Martin 2011) to remove adapter sequences. Trimmed reads with less than 30 bases were discarded. The remaining trimmed reads were aligned to the reference genome (NC_000913.3) using BWA (Li and Durbin 2010) with the -q 30 flag that trims bases from the ends of reads with qualities less than Phred score 30. SAMTOOLS version 1.3 (Li et al. 2009) was then used to generate the pileup file from the sam files generated by BWA. Finally, the list of single nucleotide polymorphisms (SNPs) and indels was compiled from the pileup file using VARSCAN version 2.3.8 (Koboldt et al. 2012). The list of mutations was filtered such that the mutations present in at least 90% of the reads, in at least one sample, were retained. Mutations present in the wildtype were filtered out.
In the RNA-seq data, ideally, the mutants should also have a mutation in yhgA as was shown in the genome sequencing done in our previous work (Mogre et al. 2014). However, the expression level of this gene is very low-only a few reads cover these genes. Since the mutated site is not covered this mutation is not called. We see the mutation in yhgA in the FusA A608E -TopA S180L second replicate though. We also see varying levels of a mutation in obgE (non-synonymous resulting in the residue change S75T) in some of the strains. This apparently mutated locus had low coverage and was mostly supported by bases towards the ends of reads. We re-analyzed the data by trimming 10 bases from both ends of the reads. After doing this, the mutation in obgE disappeared from most strains (Fig.  S2B). Thus we think that the mutation is an artifact of sequencing. We confirmed that the mutation is absent by Sanger sequencing (Fig. S2C).

Measuring variability in expression of genes within an operon
As an additional check of our RNA-seq data, we determined whether genes within operons were similarly expressed. We obtained the list of operons from the DOOR 2 database of prokaryotic operons (Mao et al. 2009). We normalized the gene read count to gene length and used this as a measure of expression of each gene. To measure the variation in expression of genes within operons, we calculated the population standard deviation (SD) of the expression values of genes within operons, individually for each operon. The formula for population standard deviation was used since a large number of operons contained only two genes (division by n instead of n-1). These operon specific measures of variation in gene expression were divided by the genome-wide variation in gene expression (i.e. the total SD across all genes for that sample, referred to as Total SD in the plots). If the within-operon SD is lower than the total SD, these values will be below 1. To check if this could arise by random chance, we randomly sampled genes equal to the number of genes in each operon. For all these sets of randomly sampled genes, we performed this calculation again. These values were also less than 1. This was because the total SD was very high, and that high variation was driven by few highly / lowly expressed genes. These genes tend to not be sampled during the random sampling because they are fewer in number. However, the median in the "Random" dataset, was much higher than that of the "Real" dataset. The difference between the two was checked using the two-sided Wilcoxon rank sum test. In all samples, the P values for comparisons between the "Real" and "Random" dataset were below 1.0 x 10 -55 .

Cell imaging and cell length estimation
Overnight grown bacterial cultures were diluted 1:100 in lysogeny broth (LB, with or without kanamycin) and grown in 96 well plates (100 µl per well) incubated at 37 °C, 200 rpm for 24 hours. Cells were pelleted by centrifugation at 4000 rpm for 3 minutes and resuspended in twice the volume of phosphate buffered saline (PBS, pH 7). These cells were then embedded under a 1% agarose pad prepared by dissolving agarose (Invitrogen) in PBS. The embedded cells were imaged using a NIKON eclipse Ti2 inverted microscope. Phase contrast images were taken using a 60× lens with oil immersion and analyzed using Oufti (Paintdakhi et al. 2016) with the E.coli_LB_subpixel.set parameters. The cell length in number of pixels was converted to microns by multiplying with a pixel to micron conversion factor (0.064). Around 725 cells of each strain were imaged from multiple experiments (n = 4) to obtain a distribution of cell lengths.

Spotting experiments to measure colony forming units (CFUs)
Overnight grown bacterial cultures were diluted 1:100 in lysogeny broth (LB, with or without kanamycin) and grown in 96 well plates (100 µl per well) incubated at 37 °C, 200 rpm for 24 hours. Serial 10-fold dilutions of these cultures in 0.75 % saline were spotted (as 5 µl spots) on LB agar plates. The plates were imaged after 14 hours of incubation at 37 °C.

Measurements of promoter activities
This assay utilizes a plasmid (pMS201) carrying a reporter gene (GFP mut2 ) downstream of the cloned promoters of genes. These plasmids were extracted from the Thermo Scientific E. coli promoter collection (Zaslaver et al. 2009). These plasmids were then introduced in the wildtype (WT with FRT site near fusA) and CyaA N600Y mutant (which also has FRT site near fusA). In this assay, GFP mut2 expression, driven by the promoters of cyaA and crp report the activities of these promoters. Simultaneous measurement of optical density at 600 nm (OD600) and fluorescence (excitation wavelength: 485 nm, emission wavelength: 520 nm) of the strains were performed in black 96-well plates, every 15 minutes, with the Tecan Infinite F200 pro plate reader, with shaking at 198 rpm and incubation at 37 C. We saw that the activities of the promoters of ⁰ cyaA and crp were higher in the mutant than the wildtype. Promoter activities were calculated as thus: Blank subtracted fluorescence values were divided by blank subtracted OD600 values. Similarly OD600 normalized fluorescence values from promoter-less strains were subtracted from these values (background subtraction, where background refers to promiscuous GFP expression from promoter-less plasmid). Promoter activity was then calculated as the temporal derivative of these fluorescence values (dF/dt). Subsequently, for each gene separately, the promoter activities in the WT and CyaA N600Y strains were scaled to fall between 0 and 1 before plotting.

Conservation analysis of RpoD
Phmmer (using the -notextw flag to facilitate sequence comparison) was used to find the best hit to the wildtype RpoD in E. coli K12 MG1655 in each of the 2074 bacterial genomes used in the analysis. Next, from each of these hits, for each residue, the number of times each residue in the hit matched the corresponding wildtype residue was counted and this count was divided by the total number of genomes analyzed to obtain a conservation score. Then we scaled the data such that it would fall between 0 and 1. This does not change the picture, only changes the scale.     Figure S5: Second site mutations, evolved in the FusA A608E background, confer growth advantage in kanamycin in the FusA P610T background as well. Growth of FusA P610T -RpoD L261Q / CpxA F218Y / TopA S180L / CyaA N600Y mutants in LB containing 8 μg/ml kanamycin (~50 % lethal concentration for the wildtype, 8-kan). Optical densities and standard deviation obtained from eight replicates is plotted. The second FusA P610T -RpoD L261Q transductant seems to have a slight growth defect which is seen in the absence of kanamycin as well (Fig. S6) and could be due to accumulation of some other mutation during the strain generation process. Figure S6: Lack of growth defect in the FusA P610T -RpoD L261Q / CpxA F218Y / TopA S180L / CyaA N600Y mutants. Growth of FusA P610T -RpoD L261Q / CpxA F218Y / TopA S180L / CyaA N600Y mutants in LB. Optical densities and standard deviations obtained from eight replicates are plotted. The domain structure of CpxA as obtained from Uniprot (http://www.uniprot.org/uniprot/P0AE82). The coloured "*"s mark the mutations in cpxA that confer kanamycin resistance.
Red "*" marks the mutation obtained in our previous work, and which we have studied in detail in this work.
Green and orange "*"s mark mutations obtained during the evolution of FusA P610T populations in 15kan, described in this work. The orange "*"s mark mutations that are present in greater than 50% of reads in at least one population, whereas the green "*"s mark mutations that are present in less than 11% of the reads. Green "*"s are SNPs. Yellow "*"s are deletions.
Plum "*" is the mutation in sbmA seen by Lazar et al. (Lázár et al. 2013).    The kanamycin resistance conferring mutation in red was obtained in our previous study and is analysed further in this study. The residues shown in orange are kanamycin resistance conferring mutations found by Lazar et al (Lázár et al. 2013). The residues shown in blue have been identified in aiding promoter escape of the RNA polymerase (Leibman and Hochschild 2007). (B) Mutated residues, and residues involved in promoter escape, are shown on the crystal structure of RpoD (PDB structure 4KMU). Colour codes of residues match that of (A). The L261Q mutation seen by us is shown in red. The L252Q mutation seen by Lazar et al. is shown in orange. The I162N mutation, seen by Lazar et al., could not be shown as the residue was absent in the structure. Correlation between the FusA A608E -RpoD L261Q mutant and the addition of fold changes of genes in the FusA A608E and RpoD L261Q mutants. In all cases the absolute of the difference of, the sum of the log2 fold changes of genes in the FusA A608E and CyaA N600Y / CpxA F218Y / RpoD L261Q single mutants, with the log2 fold changes of FusA A608E -CyaA N600Y / CpxA F218Y / RpoD L261Q double mutants is significantly greater than zero (Wilcoxon signed rank test P value < 2.2 * 10 -16 ).