Genome-wide transcriptional responses of two metal-tolerant symbiotic Mesorhizobium isolates to Zinc and Cadmium exposure

Background Mesorhizobium metallidurans STM 2683T and Mesorhizobium sp. strain STM 4661 were isolated from nodules of the metallicolous legume Anthyllis vulneraria from distant mining spoils. They tolerate unusually high Zinc and Cadmium concentrations as compared to other mesorhizobia. This work aims to study the gene expression profiles associated with Zinc or Cadmium exposure and to identify genes involved in metal tolerance in these two metallicolous Mesorhizobium strains of interest for mine phytostabilization purposes. Results The draft genomes of the two Mezorhizobium strains were sequenced and used to map RNAseq data obtained after Zinc or Cadmium stresses. Comparative genomics and transcriptomics allowed the rapid discovery of metal-specific or/and strain-specific genes. Respectively 1.05% (72/6,844) and 0.97% (68/6,994) predicted Coding DNA Sequences (CDS) for STM 2683 and STM 4661 were significantly differentially expressed upon metal exposure. Among these, a significant number of CDS involved in transport (13/72 and 13/68 for STM 2683 and STM 4661, respectively) and sequestration (15/72 and 16/68 for STM 2683 and STM 4661, respectively) were identified. Thirteen CDS presented homologs in both strains and were differentially regulated by Zinc and/or Cadmium. For instance, several PIB-type ATPases and genes likely to participate in metal sequestration were identified. Among the conserved CDS that showed differential regulation in the two isolates, we also found znuABC homologs encoding for a high affinity ABC-type Zinc import system probably involved in Zinc homeostasis. Additionally, global analyses suggested that both metals also repressed significantly the translational machinery. Conclusions The comparative RNAseq-based approach revealed a relatively low number of genes significantly regulated in the two Mesorhizobium strains. Very few of them were involved in the non-specific metal response, indicating that the approach was well suited for identifying genes that specifically respond to Zinc and Cadmium. Among significantly up-regulated genes, several encode metal efflux and sequestration systems which can be considered as the most widely represented mechanisms of rhizobial metal tolerance. Downstream functional studies will increase successful phytostabilization strategies by selecting appropriate metallicolous rhizobial partners.


Background
Metal extraction activities generate large amounts of contaminated materials. Mining spoils are major sources of pollution as they disperse into the surrounding environment through aerial or water erosion and can be found several kilometers away from their original site of deposit. Such sources of pollution pose significant risks to public health [1] and to ecosystem dynamics due to reduced biodiversity resulting from low plant coverage caused by soil toxicity [2]. While remediation through chemical extraction destructs soils and is more expensive, phytoremediation strategies appear as more environmentfriendly and constitute long term solutions to reduce metal toxicity in polluted sites. In highly polluted soils like mining soils, where the removal of metals by phytoextraction using hyperaccumulator plants is not efficient due to the slowness of the process, the best adapted method is phytostabilization [3]. This consists in limiting the dissemination of toxic metals by using metallicolous plants, i.e. metal tolerant plants, to establish a persistent plant cover and prevent pollution through erosion. Such a phytostabilization approach is nevertheless challenging since it requires metallicolous plants able to grow in soils where nutrients are most often dramatically scarce. Legume/rhizobia symbioses which transform atmospheric dinitrogen into organic nitrogen is of ecological interest here as it can improve natural soil fertility and thereby allow the colonization of other plant species and the installation of a plant cover [4,5]. However, survival and proliferation of organisms on metal contaminated sites depends on their capacity to tolerate high metal concentrations and requires the acquisition of resistance mechanisms.
A recent study conducted in a former mining area contaminated by Zinc (Zn), Cadmium (Cd) and Lead (Pb), named "les Malines Mining District" in the South of France, allowed the description of 116 plant species, some of which may be used in phytoremediation projects [6]. Among the plant species listed, a leguminous plant, Anthyllis vulneraria, is of particular interest as it can enter in symbiosis with rhizobia. A metal tolerant symbiotic partner of A. vulneraria was recently characterized as a new species of Mesorhizobium and named M. metallidurans [7]. Additionally, the Anthyllis/Mesorhizobium symbiosis has been identified in the mining site of Eylie in the Pyrenees (France). Eylie's mine soil presents a geochemical background similar to the Avinières mine where M. metallidurans was found. It is highly polluted by Zn, Cd and Pb (14,300 mg kg -1 , 23 mg kg -1 and 4,253 mg kg -1 , respectively) but displays higher organic carbon content [5] than the Avinières soil. Mesorhizobium strains isolated from Anthyllis are (i) capable of entering in symbiosis with metal-tolerant ecotypes of A. vulneraria and (ii) tolerant to several metals including Zn and Cd ; therefore, they can be used in future phytostabilization strategies and constitute a good rhizobial model for studying the mechanisms they have evolved to tolerate high metal concentrations.
Prokaryotes have developed several complex molecular mechanisms to deal with ionic homeostasis or metal toxics [8][9][10]. The most and best studied metal resistant model bacterium is Cupriavidus metallidurans CH34, which harbors in its genome an exceptionally high number of metal resistance mechanisms [11]. The major mechanisms bacteria use to counteract metal toxicity consist in (i) permeability barriers, (ii) effluxes via P IBtype ATPases, Resistant Nodulation cell-Division proteins (RNDs), Cation Diffusion Facilitor proteins (CDFs) or Major Facilitor Superfamily systems (MFSs), (iii) intracellular and/or extracellular sequestration and (iv) transformation of metals into a less toxic form by enzymatic detoxification [8]. Efflux systems represent the largest category of metal-resistance systems [9]. In the model bacterium C. metallidurans, these mechanisms are in most cases not metal specific and under the control of a complex regulatory network involving several clusters of genes and functions [12].
In order to identify genes or functions in an organism of interest, genome sequences can be produced and analyzed using comparative genomics tools. Although genomic sequence data allow for the listing of putative gene contents in an organism, they do not provide any functional evidence. To go one step further, genome-wide transcriptional analyses allow for rapid identification of genes or loci that are transcriptionally active and associated with a particular physiological state. Such transcriptomic studies can now be available for any organism at a relatively low cost using high-throughput sequencing techniques in RNA sequencing (RNAseq) experiments.
In this paper, the main goals are to study gene expression profiles associated with Zn or Cd exposure and to identify genes involved in Cd and Zn tolerance in two symbiotic Mesorhizobium isolates using comparative genomics and transcriptomics. Two strains isolated from mining spoils of distinct geographical origins are compared: strain STM 2683 was isolated from the Avinières mine in the Cévennes area [7] and strain STM 4661 was isolated from the Eylie mine in the Pyrénées Ariégeoises (France). The draft genome sequences were produced using mate-pair 454 pyrosequencing and analyzed using an automated microbial genome annotation pipeline. RNAseq analyses were performed on bacteria exposed for half a generation time to mild Zn or Cd stresses which affect growth negatively but only to a low extent. Using mild exposure allows for the detection of genes and mechanisms specifically induced in the response to metals while avoiding all stresses related to complex responses induced when submitting bacteria to high growth-inhibiting or lethal metal concentrations.

Results and discussion
Effect of Zn and Cd on Mesorhizobium metallidurans STM 2683 T and Mesorhizobium sp. STM 4661 growth The effect of Zn or Cd addition to the TY growth medium was tested for the two metal-tolerant Mesorhizobium strains (STM 2683 and STM 4661) and for the metalsensitive strains M. tianshanense (ORS 2740 T ) and Mesorhizobium sp. STM 2682 (Figure 1). Strains STM 2683 and STM 4661 were able to grow in media containing higher concentrations of Zn or Cd when compared to sensitive strains (MIC = 3 mM Zn and MIC > 0.5 mM Cd). In TY medium, the closely related Mesorhizobium strain STM 2682 or M. tianshanense were not able to grow when Zn or Cd concentrations were raised to 1 mM or 0.2 mM, respectively ( Figure 1). Aiming at characterizing genes responding to Zn or Cd exposure in our mesorhizobial strains, three RNA libraries previously depleted in rRNA for both sequenced strains were prepared. The experimental design included exposure to Cd or Zn at concentrations (0.025 and 0.5 mM, respectively) that were high enough to affect mesorhizobia growth without being lethal ( Figure 1). Such conditions were chosen in order to ascertain bioavailability of the metals and to avoid general and unspecific stress responses which would lead in the deregulation of a large number of genes that are not directly linked to the metal response. Furthermore, we used a short time exposure (approximately half a generation time) to these metals in order to detect genes and functions involved in the early response to the metals studied.

Genome sequencing, assembly and automated annotation
In order to identify the genetic bases of our two metal-tolerant mesorhizobia and to use them as references to map sequence reads from the transcriptome data sets, the draft genomes of the two studied strains were first sequenced using half a 454 run (Titanium kit) with paired-end reads from an 8 kb genomic bank for each isolate. After the assembly, the comparison with M. huakuii MAFF 303099 genome (chromosome accession number: BA000012, plasmid pMLa: BA000013 and pMLb: AP003017) led for both genomes to the identification and organization of 1 chromosomal (~6 Mb) and 1 plasmidic (~250 Kb) scaffolds. These genomes were integrated into the MicroScope platform [13] to perform automatic and expert annotation of the genes, providing 6,908 genomic objects for STM 2683 and 7,065 for STM 4661. Genome sequences are also available from European Nucleotide Archive under accession numbers CAUM01000001-CAUM01000191 for strain STM 2683 and CAAF010000001-CAAF010000089 for strain STM 4661. The general genome characteristics for the two analyzed strains are listed in Table 1. As expected, CDS located on the plasmids were less conserved. We found that only 12% and 15% of CDS located on plasmids of STM 2683 and STM 4661 respectively possessed at least one homolog in the other strain. Despite the lower conservation observed for genes located on plasmids, we did not identify an over-representation of operons or genes that were associated with metals. Such results suggest that plasmids harbor accessory genes which are not specialized in metal adaptation in the strains STM 2683 and STM 4661 contrary to the model C. metallidurans CH34 where most metal-responsive genes are plasmidborne [14].
In order to estimate genomic divergence between the two strains, we explored the proportions of homologous CDS as a function of the identity percentage. Homologous genes present in the two strains were retrieved with the comparative genomics phyloprofile tool in MicroScope. Using an alignment threshold of 80% identity over at least 80% of the query sequence and target size, approximately 4,200 homologous genes, present in both strains, were identified (4,183 CDS from STM 2683 presented at least one homolog in STM 4661 and 4,186 CDS from STM 4661 presented a minimum of one homolog in STM 2683). These results show that approximately two thirds of the genomes are conserved in the two sequenced metal-tolerant strains (data not shown).
A comparative genomic study was performed with the seven available sequenced Mesorhizobium genomes; : it shows that among annotated CDS, putative genes that encode transport and sequestration proteins were over-represented in the three metallicolous mesorhizobial genomes (data not shown).

RNAseq analysis, mapping and statistical analyses using DESeq
The universal bacterial rRNA subtraction kit (Microbe Express, Ambion) allowed for the removal of most rRNAs. Qualitative analyses of RNA samples prior and after depletion indicated that total RNA was of high quality (RIN values ranging from 7 to 10) and that most 16S and 23S rRNA peaks were removed after depletion (Additional file 1). The mRNA samples were sequenced in 36 bpcycles using the illumina HiSeq2000 (San diego, CA) with SBS technology. A lane of a FlowCell was used per sample. Image analyses and basecalling were conducted using the HiSeq Control Software (HCS 1.1.37.8) and RTA component (RTA 1.7.48). The RNAseq data are available from NCBI GEO datasets under the accession number GSE45693. The number of reads which passed the initial quality filter varied from 68.8 to ca 81.6 million and the quality analyses of the sequencing runs showed that the sequencing step was successful (data not shown). The general characteristics of the RNAseq data, listed in Table 2, showed homogeneity between treatments and strains. Notably, we found that more than 97% of the sequence reads were mapped at least once onto genomic objects identified via the automated annotation Microscope pipeline. If we consider that total bacterial RNA is composed of 95% rRNA, we found that the percentage of rRNA depletion was high (ranging from 92 to 96%) confirming that the kit used was well suited for our Mesorhizobium strains.
After removal of all rRNA reads, more than 30 million sequence reads for each treatment mapped uniquely onto genomic objects and were used to estimate gene expression. Interestingly, sequence reads were successfully assigned to all predicted genomic objects, suggesting that all CDS were transcribed and that sequencing was deep enough to cover the full transcriptomes. All read mapping data for STM 2683 and STM 4661 are accessible from the Microscope interface [15]. Using the MicroScope automated RNASeq pipeline, mapping results were then translated into raw read counts for each gene and processed through the DESeq statistical package [16] to normalize and test for differential expression between conditions (see Methods). Basically, DESeq estimates library size factors to normalize gene read counts between samples by assuming that a majority of genes have comparable expression levels in all samples. It then estimates gene expression dispersion due to biological and technical variations within each condition, and models it using negative binomial distributions. Genes for which expression levels significantly differ from the estimated dispersion are then called as differently expressed, and DESeq provides p-values adjusted for multiple testing with the Benjamini-Hochberg procedure [17] to control the false discovery rate (FDR). Although each sample (composed of 6 pooled biological replicates) was sequenced once, DESeq was able to estimate gene expression dispersion by assuming that This method overestimates dispersion relatively to when sample replicates are available, resulting in more conservative differential expression calls and slightly lower sensitivity [16]. The data produced using DESeq for all predicted CDS are accessible from the Microscope interface [15]. Possible ways to estimate the quality and reproducibility of the RNAseq data produced consist in comparing the distributions of (i) the read number per CDS prior to normalization and (ii) the differential expressions of CDS in all comparisons after normalization. Because we observed a wide magnitude in the read count numbers per CDS, a logarithmic transformation of these data was performed. A logarithmic base-2 transformation of fold changes was used to represent expression differences. The Additional file 2 reports the descriptive statistics obtained from our quantitative RNAseq datasets. The variances obtained in expression statistical analyses (Additional file 2) were low (ranging from 0.058 to 0.075), indicating that most predicted genes presented similar expression levels across treatments. Box plots can be used to compare the general characteristics of large datasets. These were produced in order to compare the distributions of our RNAseq data between treatments but also between strains ( Figure 2). The median read number per CDS was above 1,000 (log value > 3). It varied from 1,432 to 1,714 for STM 2683 and from 1,180 to 1,342 for STM 4661 (Additional file 2). The lower median read number per CDS obtained from the RNAseq data of STM 4661 can be attributed to the larger number of reads mapped on ribosomal RNAs as well as the higher CDS number that was predicted in this isolate. Nevertheless, despite this shift in the median values, comparable distributions were found as indicated by box sizes and box plot whiskers ( Figure 2A). The box plots representing the distributions of expression differences (corresponding to log2 fold changes) for all possible comparisons and for the two isolates were also very similar ( Figure 2B). In this case, the interquartile values obtained were lower than 0.3 for all comparisons, which is equivalent to a negative or a positive fold change of 1.11 or less for 50% of the annotated CDS. Furthermore, whiskers sizes (+/− 1.5 of the interquartile values representing more than 95% of the data in our distributions) were similar between comparisons and isolates and close to a log2 fold change of +/− 0.5, showing that for more than 95% of the predicted CDS, fold changes under 1.41 were observed. These results indicate that the treatments did not have a profound effect at the genomewide expression level and that the majority of CDS were not affected by the treatments. Histograms showing the distributions of the log10 read counts for each treatment and the log2 fold changes for all comparisons are presented in Additional files 3 and 4, respectively. Interestingly, they were very similar for all treatments and therefore were not strain-specific showing that the transcriptome sequencing step was reproducible at the quantitative level.
The box plot and histogram graphical representations showed that the distributions of the data were * Numbers of reads per sample obtained after mapping analyses are indicated together with their proportions relative to the total read numbers in parentheses. ** Estimated rRNA substitution rate, considering that 95% of the RNA was initially composed of rRNA. ***Only reads that uniquely mapped, and did not map rRNAs.
quantitatively comparable between treatments and strains. However, these representations do not allow to test if read numbers affect differential expression variance, especially for genes that are lowly transcribed. In order to visualize and estimate data dispersion according to read count levels, scatter plots and MA plots were produced for each Metal vs Control comparison (Additional files 5 and 6). These representations show that the data were highly centered in all comparisons. Correlation factors (R 2 values) were calculated for CDS classified by their mean read numbers for all scatter plots (Additional file 5). As expected, the correlation factor generally increased with the read count number per CDS, even though a rather high correlation factor (R 2 > 0.73) for CDSs with a mean read number under 100 was found. Altogether, our data show that the RNAseq approach produced high quality and quantitative data, even for lowly expressed CDS.

Quantitative PCR on selected genes
The reliability of our RNAseq data can be assessed by comparing the differential expression of a set of regulated genes using an alternative approach such as quantitative real time PCR (qPCR). The relative expression of five genes present in both studied microbial genomes was estimated using qPCR. Three genes presenting significant differential expression values after DESeq analyses in our RNAseq data (genes encoding for a metal-translocating P IB -type ATPase: MESS2v1_740030 / MESS4v1_360013, a periplasmic binding protein of the ABC-type transporter: MESS2v1_300037 / MESS4v1_520016, and a putative Signal peptidase II: MESS2v1_740019 / MESS4v1_360023 in STM 2683 / STM 4661, respectively) were chosen and analyzed using two reference genes (recA and glnA) whose expression was unaltered by metal treatment. The logarithm (base 2) ratios (Zn or Cd treatment / control treatment) were assessed in the two strains using the same pooled extracted RNA samples that were used for the RNAseq. The relative expression values obtained by qPCR for each gene and for the two treatments were calculated using the mean values obtained using the two reference gene expression values as standards. These relative expression levels were compared to those obtained with the RNAseq approach (Additional file 7). We found that the genes identified as significantly regulated using the genome-wide approach were also differentially expressed by qPCR for both isolates (Additional file 7A, B, D and E). Despite a higher up-regulation observed by qPCR for MESS2v1_740019 upon Cd treatment on STM 2683, a correlation factor (R 2 ) above 0.9 was found for both isolates when plotting the log2 (Metal / control) values obtained for the five selected genes using the genomewide RNAseq approach against the qPCR data (Additional file 7C and F). Differences in the amplitude of the deregulation between the two techniques could be attributed to the PCR step used to amplify libraries prior to Illumina sequencing or to differences resulting from the mRNA purification procedure. Nevertheless, the high correlation between the two approaches shows that RNAseq data can be used to assess the global gene expression levels in our Mesorhizobium strains and more particularly to identify the genes specifically induced or repressed by Zn or Cd.
To estimate the biological variance of our data, the relative expression ratios of the same target genes were also studied by qPCR in STM 2683, using total RNA isolated from independent biological replicates as starting materials, and the ratios were compared to the results obtained using pooled RNA samples (Additional file 8). Again, a high correlation factor (above 0.99) was obtained between replicates (data not shown). We found that the relative expression values obtained from the pooled RNA samples and the biological replicates were very similar with both metal treatments and the standard deviations between technical replicates were comparable to those obtained in biological replicates (Additional file 8A and B) and were highly correlated (R 2 > 0.98) (Additional file 8C). Altogether, the high correlation factors we obtained show that the RNAseq data produced in the present study contain reliable relative expression values. This is in agreement and in line with the RNAseq data analyses which showed high correlations of the read count numbers between treatments at the genome-wide level. Despite its relatively high cost, the RNAseq approach, which is based on transcript sequencing, represents a choice method for quantitative transcript measurements: the generated sequences are proofs of the corresponding gene transcriptional activities as compared to other indirect global approaches based on hybridization or qPCR which allow for the quantification of transcript levels for a limited set of targets.

Functional composition of Mesorhizobium transcriptomes and alteration upon metal exposure
The functional COG classifications of all CDS in the two Mesorhizobium genomes under study were obtained automatically using COGNiTOR [18] (Table 3).
To estimate absolute RNA contents in Mesorhizobium transcriptomes in relation to functional classes, the sums of reads for the CDS assigned to each COG functional category were calculated for each treatment. The pie charts representing the proportions of each COG class (B-W) in the transcriptomes of Mesorhizobium were produced using the means of read numbers for all three treatments ( Figure 3A and 3C for STM 2683 and STM 4661, respectively). For all the 6 transcriptomes obtained by RNAseq (2 isolates and 3 treatments), the median read number per COG class was calculated to distinguish the most represented classes at the quantitative level in the transcriptomes ( Figure 3B and 3D for STM 2683 and STM 4661, respectively). The most represented functional categories at the transcriptional level in both isolates belonged to metabolism processes, principally COG classes E, G, P and C which refer to the metabolism and transport of amino acids, carbohydrates, inorganic ions and energy production and conversion, respectively, the poorly characterized COG class R, and class J which refers to translation and ribosome structure and biogenesis. Our quantitative analysis highlights the relative content of each functional class in the transcriptomes likely to be representative of Mesorhizobium growth in the medium under study. However, as expected, the number of transcripts per COG class was correlated to the number of CDS per class (R 2 of 0.8 data not shown) and is thus reflecting the relative genomic content. In order to estimate the relative transcriptional activity of each functional category, the sums of reads per COG class was normalized to their respective CDS numbers ( Figure 4). Such normalization shows that the COG classes belonging to (i) cellular processes and signaling (classes M, O, T, U, N and D) (ii) information and storage processing (classes J and L) and (iii) metabolism (class C) could be considered as transcriptionally over-active as compared to the other COG classes in both isolates ( Figure 4 and Table 3). These results are not surprising as the cells were harvested during the exponential growth phase which requires the involvement of many genes whose function is associated with (i) basic cellular processes such as cell cycle control, cell division, chromosome partitioning, cell wall/ membrane/ envelope biogenesis, signal transduction mechanisms, intracellular trafficking and (ii) information storage and processing such as replication and translation, all of which require high energy production and conversion.
Interestingly, slight differences were noted for the COG class J that refers to translation, ribosomal structure and biogenesis. Indeed, in both isolates, the two metal treatments triggered a decrease of the transcription levels for genes assigned to this particular COG class ( Figure 3 and 4). Box plot representations of the log2 fold changes between metal and control treatment data as classified by COG functional categories (Additional file 9) confirmed this observation, suggesting that metal treatment affects translation and ribosome structure and biogenesis. These results are in line with observations made by Pereira and colleagues (2006) who found that metal exposure could affect protein levels as shown for Rhizobium leguminosarum bv. viciae isolates [19]. However, despite a relatively high tolerance of our Mesorhizobium isolates as compared to other mesorhizobia [7], our transcriptomic data suggest that a global repression of the translation machinery occurs upon Cadmium or Zinc exposure; according to Pereira's hypothesis, such repression should categorize our isolates as metal sensitive. Additional experiments are required to verify whether our isolates are affected in protein synthesis and compare their levels to other mesorhizobia, including strains more tolerant than those presently studied in order to validate this hypothesis and show that it is applicable to the Mesorhizobium genus. Interestingly, previous studies on E. coli show that upon Cd exposure, protein biosynthesis machinery was stopped [20].
In order to assess the effect of metals on general functional processes in the isolates, we compared the functional distribution of the most significantly deregulated genes (padj ≤ 0.1) to which a COG class was assigned to full gene set distribution ( Figure 5). The percentage of deregulated genes per category gives a rapid mean to estimate which COG classes are the most affected by metal treatment. Interestingly, classes O and P appeared as over-represented when compared to the full gene set for both metals and both isolates. Class O refers to post-translational modification, protein turnover, chaperones, which supports the hypothesis that protein biosynthesis and stability are affected by Zn and Cd exposure. Class P corresponds to inorganic ion transport and metabolism, which indicates that Zn or Cd divalent cations are processed by bacteria as soon as they are sensed, most probably through export systems [9].

Identification of genes regulated by metals in Mesorhizobium strains STM 2683 and STM 4661
Transcriptional profiles obtained using the genome-wide RNAseq approach was performed on pooled samples. The high reproducibility between treatments, as exemplified by the high correlation factors of read counts and relative expression values of the RNAseq datasets allowed us to obtain adjusted p-values (FDR) with DESeq. Putative CDS presenting FDR values < 0.1 in a "Metal" vs "Control" experiment are listed in Tables 4 and 5 for STM 2683 and STM 4661, respectively. A total of 72/6,844 and 68/6,994 putative CDS (or 1.05% and 0.97% CDS for STM 2683 and STM 4661, respectively) were found to be significantly differentially expressed by one or both metals for STM 2683 and STM 4661, respectively. These putative differentially expressed CDS are classified into six categories: transport, sequestration, regulation, oxidoreduction, others and unclassified. Thirteen out of seventy-two and 13/68 putative CDS (18% for STM 2683 and 19.1% for STM 4661) were unclassified for STM 2683 and 4661, respectively. As expected, the category "unclassified" is the most represented among significant regulated putative CDS. Few CDS are classified into the category "others" which can be assimilated to the non-specific metal response, suggesting that the RNAseq approach in our conditions is reliable and allowed us to identify the genes involved in the specific metal response. "Transport" and "sequestration" are the most represented categories after the "unclassified" category among significantly regulated putative CDS, which is in agreement with the fact that transport and sequestration proteins were found over-represented in the three available metallicolous mesorhizobial genomes. Among the putative CDS significantly regulated in STM 2683, 39 ( (Tables 4 and 5). Among the 1.05% and 0.97% putative CDS that are significantly and differentially expressed for STM 2683 and STM 4661, respectively, 13 putative CDS presenting homologs for both isolates were identified as significantly regulated under Zn and/or Cd treatments. Among these conserved and significantly regulated putative CDS, some  Table 3 for description) in the Mesorhizobium transcriptomes. Histograms representing the quantitative composition of the transcriptomes (B and D) organized by the biological Process and Class ID (bottom of histograms) show the total amounts of reads (indicative of transcript abundance) in each class and for each treatment (Zinc, Control or Cadmium). The median read numbers per class and per treatment are indicated as the last bars on the histograms which were used to draw the dotted line above which the defined classes are the most represented for both organisms and for all treatments (these class IDs and their proportions are indicated in white on the pie charts).
had homologies to genes encoding proteins whose function could be associated to specific metal responses like metal sequestration or metal transport.

Conserved and differentially regulated genes possibly involved in sequestration
Metal sequestration can be used by microorganisms to reduce metal bioavailability and thereby avoid metal toxicity. Metal sequestration includes cell wall components, periplasmic or cytoplasmic binding proteins involved in the precipitation of metals by phosphates or in their reduction as sulfides, and intracellular binding proteins containing CPX-Zn finger domains.
Among the 13 putative CDS present and regulated by metal(s) in both isolates, several could be involved in extracytoplasmic sequestration. Indeed, we identified a putative opgC gene encoding a membrane protein required for succinylation of periplasmic glucans in M. loti [21]. It was slightly but significantly induced upon Cd exposure (> 1.6 fold) in both Mesorhizobium strains (MESS2v1_310036 in STM 2683 or MESS4v1_520093 in STM 4661) but was not significantly regulated by Zn (Tables 4 and 5). The addition of succinyl groups confers global negative charges to periplasmic cyclic glucans and may thus limit the diffusion of toxic cations into the cytoplasm through immobilization in the bacterial cell membrane. Amongst the conserved and regulated CDS that may be involved in metal precipitation, we also identified a putative signal peptidase II, also referred to as aspartic peptidase, and MEROPS family A08 (MESS2v1_740019 in STM 2683 and MESS4v1_360023 in STM 4661) significantly up-regulated upon Zn exposure (> 3 fold) in both isolates. It presents similarities to the polyprotein peptidase PbrC of C. metallidurans strain CH34 that could participate in the precipitation of Pb 2+ via the generated phosphate ions thus limiting its re-entry into the cell cytoplasm after extrusion via the P IB -type ATPase PbrA [22]. Interestingly, a putative membraneassociated phospholipid phosphatase (MESS2v1_740015 in STM 2683) located in the vicinity of the putative signal peptidase II, aspartic peptidase, MEROPS family A08 was also significantly up-regulated by Zn in STM 2683 (> 2 fold). Precipitation of metals such as Cd by surface polysaccharides, phosphates, sulfides or others to form insoluble salts that reduce metal bioavailability has been observed in several microorganisms [23][24][25]. Additionally, RNAseq data allowed us to identify a few proteins with a thioredoxin domain significantly up-regulated by Cd in STM 2683 (MESS2v1_1660012) and up-regulated by both metals in STM 4661 (MESS4v1_360050). Thioredoxin is a general protein disulfide reductase believed to serve as a cellular antioxidant by reducing protein disulfide bonds produced by various oxidants and also interact with other proteins to form functional protein complexes [26]. A recent study has demonstrated the role of thioredoxin in the reduction of U (VI) to U(IV) and Cr(VI) to Cr(IV) in Desulfovibrio desulfuricans G20 [27]. Among our RNAseq data, we identified other CDS significantly up-regulated in the two Mesorhizobium strains under Zn and/or Cd exposure (7 CDS identified in STM 2683 and 8 in STM 4661) that may be involved in the precipitation of metals by phosphates or their reduction as sulfides (Tables 4 and 5).
Among the putative CDS possibly involved in metal sequestration, we also identified several conserved hypothetical proteins that contained two cytoplasmic CPX-Zn finger domains (MESS2v1_730267 in STM 2683 and MESS4v1_360040 in STM 4661) that were significantly up-regulated (> 2 fold) upon Zn exposure in both Mesorhizobium strains. CPX-Zn finger domains are known to coordinate Zn ions with a combination of Figure 5 Impact of metals on the major functions. Histograms representing the percentages of CDS classified by functional COG category for STM 2683 (A) and STM 4661 (B) in the total set of CDS (black bars), in the Zinc or Cadmium treatments when compared to the control treatment (respectively light and dark grey bars). For a defined COG class, grey bars above the black ones indicate that upon metal treatment, the corresponding COG class contains a higher proportion of CDS that display differential regulation.  cysteine and histidine residues. Interestingly, another homolog was also found in STM 2683 (MESS2v1_ p180003), which was significantly induced by both metals (> 4 fold by Zn and > 3 fold by Cd). Several studies suggest a role of cytoplasmic proteins such as Zn-finger proteins or cysteine residues in metal binding. On the cellular speciation of Cd ions in E. meliloti, if extracellular sequestration was favored, a significant amount of intracellular metal was however measured [28] which most probably results from chelation. The transcriptomic analyses of E. coli in response to a Zn stress showed that adding excess external Zn induced the expression of many genes that are organized in the regulon for cysteine biosynthesis, implying that Zn-binding proteins containing cysteine residues are one of the mechanisms that chelate Zn [29]. Finally, the metallothionein SmtA and several SmtA-like proteins have been described to sequestrate and detoxify Zn and Cd through Zn-fingers in Synechococcus PCC7942 , in P. aeruginosa and P. putida [30].

Conserved and differentially regulated genes possibly involved in Transport
Metal transport can be used by microorganisms to limit metal entry and thereby limit negative effects on cell components. Metal transport involves efflux proteins such as Resistant Nodulation cell-Division proteins (RND), P IB -type ATPases, Cation Diffusion Facilitor proteins (CDF) or Major Facilitor Superfamily systems (MFS) [9], multi-drug transporter proteins and proteins involved in metal homeostasis.

Identification of candidate genes possibly involved in transport by CDF efflux systems
Among the putative CDS probably involved in metal efflux, we identified a putative CDF-type transporter that presented homologies with ZitB, a Zn-specific transporter identified in E. coli [31] and CzcD involved notably in Cd and Zn resistance in C. metallidurans [32] that was significantly up-regulated (> 2-fold) after Zn treatment in STM 4661 (MESS4v1_360039). CDF transport systems function as cation/proton antiporters and are driven by concentration, chemiosmotic gradient or potassium gradient [33]. Surprisingly, in STM 2683, the homolog was fissed into two separate CDS (MESS2v1_730268 and MESS2v1_730269) that were similarly regulated. Whether these two CDS were still functional in STM 2683 is not known. We also found in STM 2683, a second CDS (MESS2v1_p180001) showing strong homologies to the 5' end of CzcD that was induced by Zn (4-fold) and Cd (> 2fold). However, this particular CDS was located in the edge of a contig, and we therefore lack its 3' end indicating that it could actually be the true homolog of MESS4v1_360039.

Identification of candidate genes possibly involved in co-transport of metals and antibiotics
Additionally, we identified two putative CDS mdtC and mdtA, encoding a permease protein and a membrane fusion protein, respectively (MESS2v1_620035-36 in STM 2683 and MESS4v1_690033-34 in STM 4661). They were significantly induced by Zn in STM 2683 but their homologs in STM 4661 were not significantly up-regulated by metals. Both genes encoded putative RND-type proteins involved in multidrug efflux. In Lactococcus lactis and E. coli, MdtA confers multiple antibiotic resistances [34]. It was also reported that in Salmonella, mdtABC contributes to Cu and Zn resistance, in addition to their role in the resistance to β-lactams, novobiocin and deoxycholate [35]. In E. coli, a microarray analysis also demonstrated that the mdtABC operon was up-regulated in response to stress caused by excess Zn, suggesting that metal ions can modulate bacterial resistance   to antibiotics [36]. Possible explanations underlying the coselection process between metals and antibiotics include co-resistance (different resistance determinants present on the same genetic element) and crossresistance (the same genetic determinant responsible for resistance to antibiotics and metals) [37].

Identification of candidate genes possibly involved in transport by P-type ATPase systems
Among the 13 conserved and regulated CDS probably involved in metal efflux, we identified two P IB -type ATPases. P IB -type ATPase transporters belong to the large superfamily of ATP-driven pumps involved in the transmembrane transport of charged substrates against their concentration gradients [38]. P IB -type ATPases are known mechanisms classically used by microorganisms to maintain constant and non-toxic levels of metals in the cytoplasm. Transcriptomic analyses of E. coli or E. meliloti under Zn and Cd stress showed the high induction of genes encoding P IB -type ATPases [36,39]. Exclusion from cells by efflux transporters like the P IB -type ATPase ZntA of Escherichia coli [40] or CadA of Staphylococcus aureus [41] are among the most studied bacterial P IB -type ATPases. We identified, in the RNAseq data of the two Mesorhizobium strains, a gene encoding a related ZntA (MESS2v1_980073 in STM 2683 and MESS4v1_610026 in STM 4661) which was significantly induced upon Zn exposure (> 25-fold in STM 2683 and > 3-fold in STM 4661). Intriguingly, the predicted ZntA was also strongly induced by Cd (>11-fold) in STM 2683, while its homolog in STM 4661 was not up-regulated, indicating that different regulatory mechanisms could operate in the two strains. In accordance with this assumption, a putative MerR-type regulator with homologies to ZntR was predicted upstream of ZntA in STM 2683, on the opposite strand. Furthermore, we found highly different numbers of putative metal responsive regulators (HmrR) predicted in their genomes. A total of 19 possible HmrR-encoding CDS were predicted in STM 2683, including a ZntR and a MerR homolog which localized close to a locus involved in Mercury transport and reduction. In contrast, only four putative HmrR-encoding genes were predicted in STM 4661, none of them being automatically annotated as ZntR. The large number of predicted HmrR-encoding genes present in STM 2683 suggests a more complex regulation network in this strain when compared with STM 4661 which could explain the different transcription profiles obtained for several metal regulated genes in the two Mesorhizobium strains. Another P IB -type ATPase . ccdA encodes a thioldisulfide oxidoreductase involved in the cytochrome c-type biogenesis [42]. A dsbD-related gene has been identified as a metal-binding motif involved in Copper tolerance in E. coli [43]. The precise role such homologs have in our Mesorhizobium strains remains to be addressed even though it is tempting to speculate on a possible involvement in metal binding via divalent metal cation-sulfide bond formation or in the functioning of the P IB -type ATPase by providing a reducing power.
Identification of candidate genes possibly involved in transport by ABC transport systems Among the differentially regulated CDS, we identified three conserved CDS that are organized into an operon (MESS2v1_300039-37 and MESS4v1_520018-16 for STM 2683 and STM 4661 respectively) encoding for an ABC-type transport system. These genes were significantly repressed after Zn exposure (> 17-and > 6-fold repression in STM 2683 and STM 4661, respectively) and induced upon Cd exposure (>5-and >9-fold induction for STM 2683 and STM 4661, respectively). These three CDS contain conserved domains related to ABCtype Mn 2+ /Zn 2+ transport systems, including respectively the multidomain ATPase component ZnuC (COG1121), the membrane permease ZnuB (COG1108) and the soluble periplasmic Zn-binding protein ZnuA (cd01018) that captures Zn and delivers it to ZnuB. These znuCBA orthologs are strongly repressed upon Zn exposure ; therefore they may be involved in Zn homeostasis in our isolates, as showed in E. coli [44]. The transcriptional regulation of this high affinity Zn-import system has been shown to be under the control of a Zur (Zn 2+ Uptake Regulator) repressor that belongs to the Fur (Ferric iron Uptake Regulator) family of transcriptional regulators [44,45]. Interestingly, orthologs of M. loti and E. meliloti Zur proteins are present in STM 2683 and STM 4661 (MESS2v1_1530029 and MESS4v1_830491, respectively) and Agrobacterium-specific conserved zur binding sites [46] were also found in the promoter region of the znuCBA putative operon in both mesorhizobial strains (data not shown). When bacteria were subjected to Zn in the medium, the znuCBA genes were strongly repressed, suggesting, together with high homologies to other bacteria, that the presently identified ATP-binding cassette transport system is also a high affinity Zn-import system involved in Zn homeostasis in M. metallidurans. Moreover, previous work showed that other divalent cations such as Fe 2+ , Mn 2+ , Cu 2+ or Cd 2+ de-repressed the expression of the znu operon in E. coli [45]. More recent studies showed that the expression of znuA and znuB increased upon Cd exposure in Salmonella sp. [47] and that ZnuAtagged proteins increased upon Cd exposure in E. coli O157:H7 [48]. This regulation could either be direct (by direct binding of divalent cations to Zur, thereby changing the conformation of the regulator) or indirect (by lowering the intracellular Zn levels), as suggested previously [45,47]. A possible explanation for znu induction by Cd could be that direct binding of divalent Cd cations to the Zur regulatory repressor protein as a substitute for Zn alters its conformation and makes it unable to bind DNA, thus directly limiting its repressive activity. Another possibility could be a deficit in intracellular Zn concentration most probably through competition between Zn and other cations for transporters [45,47]. In line with this later observation, we found that both Zn and Cd induced several efflux systems in our isolate (Tables 4 and 5) that may participate in the unspecific export of divalent cations, such as the P IB -type ATPase CadA (encoded by MESS2v1_740030 in STM 2683 and MESS4v1_360013 in STM 4661), which is significantly induced by both metals and appears important for Zn adaptation in STM 2683 (Maynaud et al., unpublished observations). To increase our knowledge about the regulation of this high-affinity Zn-import system, we used a promoter-PROBE reporter system [49] and analyzed the transcriptional activity of the znuC (MESS2v1_300037) in the presence of various metal concentrations and mixes in the wild type STM 2683 and in the cadA deleted mutant. We showed that the mutant was affected in Zn tolerance as previously described (Maynaud et al., unpublished observation). We tested whether a reduction of the efflux components in STM 2683 through cadA mutation could alter the regulation of znuC upon exposure to Zn/Cd or to various mixes of the two metals ( Figure 6). Increasing the Zn concentration in the medium reduced significantly the relative expression of znuC in the wild-type bacterium, reaching non detectable levels when 0.4 mM Zn was used, while no transcriptional activity was detected in the cadA mutant ( Figure 6A). Since no Cd was present in the control TY medium, the fact that no induction was observed in the cadA mutant, while being relatively highly induced in the wild type, indicates that intracellular Zn concentrations were high enough in the mutant to repress znuC transcription, and this may result from its lower efflux capacity. When Cd was added to the growth medium, we found that znuC expression increased in the wild-type bacterium in the presence of as low as 0.00625 mM, but higher Cd concentrations did not further increase its expression levels and even tended to decrease them when Cd concentration were above 0.2 mM ( Figure 6B). Surprisingly, znuC expression required higher Cd concentrations in the mutant, reaching a maximum at 0.1 mM Cd and then decreasing again. These results support a more indirect effect of Cd ions which could reduce intracellular Zn levels, leading to a higher stimulation of the high-affinity Znimport system. To test whether Cd competed directly with the Zur regulator, the effect of various Cd/Zn mixes on znuC transcriptional activity was studied ( Figure 6C). No transcriptional activity was observed in the cadA deletion mutant suggesting that intracellular Zn concentrations were high enough in that strain to repress znuC transcription and out-competed the Cd level which, furthermore, was probably extruded or sequestered via other mechanisms. Yet, interestingly we found that the molar ratio of Zn and Cd ions in the medium affected znuC transcription in the wild-type bacterium. Indeed, despite a metal-dose dependent transcriptional activity of the efflux pump CadA (Maynaud et al., unpublished observation), increasing Zn/Cd molar ratios from 0.5 to 2 negatively affected the transcriptional activity of znuC in STM 2683, regardless of the total metal concentration used. These results suggest that the two divalent cations compete for binding to the Zur regulator which consequently alters the regulation of znuC. Our znuC transcriptional data indicate that Cd generates both direct and indirect effects on znuCBA regulation in STM 2683, but additional experiments are required to decipher this complex regulation and the effects of Cd on Zn homeostasis.

Conclusions
In the present study, we used comparative genomic and transcriptomic analyses to study the gene expression profiles associated to Zn and Cd exposure and to identify putative genes involved in metal tolerance in two metallicolous Mesorhizobium strains isolated from Anthyllis nodules from two distinct mining sites in France. A RNAseq-based approach allowed for the rapid discovery of homologous genes that responded specifically to metals in both strains.
Analyses of the genes identified as significantly regulated suggest that transport and sequestration are the most important mechanisms underlying the metal response in metallicolous mesorhizobial strains. We identified proteins involved in metal homeostasis (ABC-type transport system ZnuABC), metal or multi-drug efflux systems which allow for intracellular detoxification (CadA, related ZntA, related CzcD/ZitB, related MdtC and MdtA) and several cell membrane components, periplasmic and cytoplasmic proteins possibly involved in metal precipitation and binding to reduce metal bioavailability. Surprisingly, the number of HmrR regulators and the inductibility of ZntA differed between our two Mesorhizobium strains. Additionally, the comparative RNAseq-based approach revealed a low number of genes significantly regulated in both strains (around 1%) and a low number of genes involved in the nonspecific metal response, indicating that the approach was fit for identifying genes that specifically respond to metals. Global RNAseq analyses suggested that Zn and Cd also repressed the translational machinery. Our study allowed to detect genes that are conserved among two metallicolous mesorhizobia isolated from distant sites, whose expression is metal-dependent and which show a potential involvement in metal tolerance. Studies on the distribution of potential genetic determinants identified in the present work in other metal tolerant and sensitive Anthyllis mesorhizobia will enable the definition of biomarkers. Such candidates may help us to select appropriate rhizobial symbionts to be used as biofertilizers to improve phytostabilization strategies and limit the impacts of mine spoil dispersal on human health and surrounding environments. Finally, functional studies on the molecular determinants identified herein and their speciation during in situ symbiotic growth will allow to evaluate their implication in the metal tolerance of our isolates and their role in symbiotic plant growth promotion on contaminated sites.

Microbiological methods
The bacterial strains and plasmids used in the present study are listed in Additional file 10 [5]. Mesorhizobium tianshanense type strain ORS 2640 T was used as a metal-sensitive control.

STM 2683 and STM 4661 genome sequencing, assembly and automated annotation
Genomic DNA from strains STM 2683 and STM 4661 were prepared from 15-ml cultures following previously described protocols [52]. The complete genome sequences of strains STM 2683 and STM 4661 were determined using 454 Titanium technology. A mate-pair genomic library of 8 kb was constructed using Roche reagents from genomic DNA fragmented with HydroShear system. Then sequencing was carried out on this library, giving reads of 350bp average length, up to a final coverage of 30× for STM 4661 and 33× for STM 2683. Assembling was done using Newbler (2.3 version), and a comparison with Mesorhizobium loti genome was performed to identify the scaffold organization. The draft genome sequences were uploaded into the MicroScope platform [13] and subjected to the automatic annotation pipeline. The genome sequences are available and can be browsed at the MaGe interface [53] and are part of the "Rhizoscope" project [54].  grown in TY medium up to mid-exponential phase (up to OD 600nm = 0.5) in 250-ml Erlenmeyer flasks at 28°C and 145 rpm. When the pre-cultures reached the midexponential phase (after ca 16 hours), 15 ml were added to an equal volume of pre-warmed TY containing either nothing (Control treatments), Zn at 1 mM (Zn treatments, 0.5 mM final concentration) or Cd at 0.05 mM (Cd treatments, 0.025 mM final concentration) and were further incubated for 5 hours, which corresponds to ca half a generation time, at 28°C at 145 rpm. After 5 hours incubation, 1/10 th volume of ice-cold stop buffer (5% phenol in ethanol) was added to each culture and directly centrifuged at 4°C, 8,000 rpm for 4 minutes. Supernatants were quickly discarded and the tubes were re-centrifuged for 1 minute so as to collect all the liquid at the bottom of the tubes which was discarded by pipetting. Cells were snap-frozen in liquid nitrogen and stored in the −80°C freezer until use.
Total RNAs were purified from each replicate, each treatment and each strain individually using the RiboPure™ kit (Ambion) following the manufacturer's recommendations. After isolation of total RNA, a DNase I treatment was performed as described in the RiboPure™ Bacterial kit protocol. Equal amounts of total RNAs from the biological replicates belonging to the same strain and the same treatment were pooled and subjected to two successive runs of ribosomal RNA subtractions using the Microbe Express™ kit (Ambion) following the manufacturer's instructions. The rRNA removal procedure consists in incubating RNA samples with magnetic beads to which complementary 16S and 23S rRNA conserved sequences are linked. After hybridization, captured 16S and 23S rRNAs are pulled to the side of tube with magnets allowing unbound mRNAs to be eluted. After two rounds of subtraction, mRNA samples were finally precipitated with ethanol and resuspended in nuclease free water. Total and messenger RNA quantities and quality were assessed by Nanodrop and Bioanalyser using RNA chips for Agilent 2100 Bioanalyser (Additional file 1).
For RNA sequencing, six cDNA libraries corresponding to the three treatments (control, Zn and Cd), 100 ng of mRNA were used according to the Illumina kit instructions. Briefly, the first step in the workflow involves fragmentation of mRNA into small pieces using divalent cations under high temperature. Then the cleaved RNA fragments were copied into first-strand cDNA using reverse transcriptase and random primers. This was followed by second-strand cDNA synthesis using DNA Polymerase I and RNase H. These cDNA fragments then went through an end-repair process, the addition of a single ' A' base, and then ligation of the adapters. These products were then purified on a gel to select a size range at 200 pb and enriched by PCR to create the final cDNA library. The efficacy of the library construction was checked in a quality control step that involved measuring the adapter-cDNA size and concentration on an Agilent DNA 1000 chip. Sequencing libraries were denatured with sodium hydroxide and diluted to 6 pM in hybridization buffer for loading onto a single lane of an Illumina HiSeq 2000 flowcell V1.5. Cluster formation, primer hybridization and single-read, 36 sequencing cycles were performed on cBot and HiSeq2000 (Illumina, San Diego, CA) respectively.
After sequencing, the sequences or reads that passed the quality filter were mapped on the annotated genomic objects identified in the genomes of STM 2683 and STM 4661 using ssaha2 software [55]. Mapped reads were converted to read counts per gene using BEDtools [56], which were then analyzed using DESeq standard protocol for conditions with no replicate [16]. DESeq generated differential expression values (fold changes) for all genomic objects and all pairwise comparisons of conditions and associated adjusted p-values (padj) controlling for the false discovery rate (FDR) to them. Raw and analyzed RNAseq data are accessible at the MicroScope web interface [15] (Additional file 2).
The distributions of (i) the raw read count numbers per CDS after log10 transformation for each treatment (Additional files 3 and 5) and (ii) the fold changes after normalization and log2 transformations for each intraspecies comparison produced by DESeq (Additional files 4 and 6) were analyzed using Microsoft Excel and XLSTAT TM visualization tools such as box plots, histograms, scatter plots or pie charts. The functional COG classification of all CDS in the two genomes was done automatically using COGNiTOR which is included in the automated annotation MicroScope pipeline. CDS with assigned COG classes were retrieved from Micro-Scope and used to estimate the general composition of the transcriptomes of the two strains as well as the overtranscribed classes. The global effects on the major functional categories Zn and Cd treatments had on the two strains were assessed by comparing the distributions of CDS fold changes classified by COG categories and by comparing COG classes to which the most significantly regulated genes (FDR < 0.1 in at least one comparison) were associated (Additional file 9).

Real time PCR on selected genes
The primers used for quantitative PCR experiments are listed in Additional file 11. For each technical replicate, 1μg of total RNA from either total pooled RNA used for RNAseq or total RNA isolated from independent biological samples prior to pooling were reverse-transcribed with 400 U of Super-Script II (Invitrogen) and random hexamer primers. qPCRs were performed on a Strategene MXP3005P system using Power SYBER green master mix (Applied Biosystems). PCR started with 10 min incubation at 95°C, followed by 40 cycles consisting of 15 seconds at 95°C and 30 seconds at 60°C. Primer specificity and the formation of primer dimers were monitored by dissociation curves. The expression levels of the metaltranslocating P IB -type ATPase encoding genes (cadA) (MESS2v1_740030 and MESS4v1_360013), the periplasmic binding protein of the ABC-type transporter genes (znuA) (MESS2v1_300037 and MESS4v1_520016) and the putative Signal peptidase II genes (MESS2v1_740019 and MESS4v1_360023) were standardized by using recA (encoding part of the DNA recombination and repair system) (MESS2v1_330003 and MESS4v1_540005) and glnA (encoding part of the glutamine synthetase I) (MESS2v1_390006 and MESS4v1_580041) as references. PCR efficiency (E) for each amplicon was calculated using the linear regression method on the log (fluorescence)per-cycle-number data using Stratagen MXpro software. All qPCRs were performed in three technical replicates using either total pooled RNA used for RNAseq (two replicates for each mesorhizobial strain) or total RNA isolated from three biological samples prior to pooling so as to estimate biological variability (three replicates on STM 2683). Data for each sample were expressed relative to the expression levels of recA or glnA by using the mathematical model described previously, which determines the relative quantification of a target gene in comparison to a reference (ref) gene between treatment and control samples [57]. The relative expression ratio (R) of a target gene is calculated based on efficiency (E) and the threshold cycle (Ct) of an unknown sample versus a control and expression in comparison to a reference gene. The average relative quantity for each gene under study was calculated and log2-transformed. Means and standard deviations of these final log2 ratios were calculated using data from three technical replicates. Mean comparison was performed with a multiple comparison of Conover-Iman, using XLSTAT software, to determine significant differences in gene expression levels between inducible genes and reference genes and to determine significant differences between the metal treatments (Cd or Zn condition) and the control metal-free treatments.

Use of promoter-gfp fusion in transcriptional analyses
The promoter region of the ABC transporter ATP-binding protein encoding gene znuC (MESS2v1_300039) was amplified using primers STM 2953-Forward and STM 2954-Reverse which generate a 500-bp PCR product. PCR products were cleaned up using illustra™ GFX™ PCR DNA and Gel Band Purification Kits (GE Healthcare). The purified PCR product was ligated into a pGEMt-easy vector (PROMEGA) following the manufacturer's recommendations to generate plasmid pGEMt-znuCp that were transferred into E. coli XL2 Blue Ultra-competent Cells. Inserts from pGEMt-znuCp were excised using the restriction enzyme HindIII and ligated into a pPROBE GT vector which allows for transcriptional fusions with the gfp gene [49], previously digested with HindIII and dephosphorylated, and the ligation mix was transferred into electrocompetent E. coli strain DH10B. Plasmid DNA was extracted using the Wizard Plus SV Minipreps DNA Purification System kit (Promega) following the manufacturer's recommendations. Plasmid containing the promoter of znuC (pGT-znuCp) was sequenced to verify sequence integrity.
The construction (pGT-znuCp) was transferred into STM 2683 electrocompetent cells as described below. Electro-competent M. metallidurans STM 2683 T cells were prepared following a method developed previously [58] but with slight modifications. Briefly, 5 ml of a fresh bacterial culture were used to inoculate 100 ml of TY medium and incubated at 28°C with shaking at 250 rpm until OD 600nm reached 0.5 to 0.7. Cells were harvested by centrifugation at 10,000 rpm for 5 minutes at 4°C. The pellet was washed with 10 ml of ice-cold 300 mM sucrose solution and centrifuged as before. The pellet was finally resuspended in 5 ml of ice-cold 300 mM sucrose containing 15% glycerol. Aliquots of 100-150 μl of electro-competent cells were snap-frozen in liquid nitrogen and stored at −80°C until use. When needed, electro-competent cells were transformed using the following procedure. Cells were thawed on ice and placed in pre-cooled 0.2 mm electroporation cuvettes for 20 minutes on ice. Five to ten μl of 30-50 ng μl -1 of plasmid DNA were added to the cuvettes and electroporation was conducted at 12.5 KVolts and 720 ohms. After pulse application, the cell suspensions were diluted with 1 ml of TY medium and incubated at 28°C with shaking for 11hours. After incubation, 100 μl of cells were plated onto TY plates containing appropriate antibiotics. Transformants were selected as single colonies and the presence of the plasmid was verified by PCR.
When needed, STM 2683 carrying the promoter-pPROBE constructs were grown on TY agar plates supplemented with antibiotics. Using sterile toothpicks, fresh colonies were transferred to sterile 8-tube strips containing 100 μl of TY. Cells were homogenized by pipetting and for each transcriptional assay, the same bacterial inocula (10 μl) were used to inoculate 150 μl TY medium supplemented or not with Zn (0.0125, 0.025, 0.05; 0.1, 0.