Nucleotide sequences of heavy metal P_1B-ATPases showing orthology to already identified heavy metal P_1B-ATPases were extracted from the Genbank sequence databases following similarity analyses performed using the TBLASTN program with the filtering removed. The sequence databases that were analysed were the ‘non-redundant’ (nr) database, the ‘high throughput genome sequences’ (htgs) database, the ‘whole-genome shotgun reads’ (wgs) database and the ‘NCBI genome’ (chromosome) database. The analyses were completed considering the sequences available in databases at the end of April 2011. The set of sequences that we obtained is extensively described in Table 1.

To deduce the protein sequences from the raw nucleotide sequences, annotation was performed mainly through similarity analyses. The query DNA sequences were compared: (i) to protein sequences translated from validated HMA encoding cDNAs in pairwise comparisons using TBLASTN; as well as (ii) to genome sequences from other species using TBLASTX. The precise positioning of the introns was performed through the recognition of the GT and AG motifs delimiting the 5’ and 3’ ends of introns, respectively. This approach proved to be more precise and reliable than using annotation softwares such as FGENESH, Eukaryotic GeneMark, EuGène or SpliceMachine.

To determine the length of the C-terminal cytoplasmic domain, we started from the analysis performed for the AhHMA4 gene [12], which indicated the position of the transmembrane helices. Using both a sequence similarity approach and the TMHMM software (http://www.cbs.dtu.dk/services/TMHMM/), we determined the position of the last transmembrane domain in all the HMA sequences. We then considered that the C-terminal loop started at the first amino acid located downstream the last transmembrane domain.

For the phylogenetic analyses, full-length amino acid sequences were aligned by CLUSTALW and imported into the Molecular Evolutionary Genetics Analysis (MEGA) package version 4 (http://www.megasoftware.net [24]). All positions containing gaps and missing data were eliminated from the dataset (Complete deletion option). Phylogenetic analyses were carried out using the Neighbor-Joining method [25], the Minimum Evolution method [26], the UPGMA method [27], the Maximum Likelihood method [28] and the Maximum Parsimony method [29], along with statistical bootstrapping procedure involving 1000 replicates.

In order to perform the phylogenetic analysis of the Zn/Cd/Pb/Co P_1B-ATPases family, it was decided to establish a complete set of sequences. The A. thaliana HMA1, HMA2, HMA3 and HMA4 protein sequences – which were deduced from the nucleotide sequences whose accession numbers were NM_119890, AY434728, AY055217 and AF412407, respectively – were used as starting queries to search for nucleotide sequences encoding similar proteins in the public sequence databases using the TBLASTN program. All the different databases gathering genome sequences were analysed. The ‘expressed sequence tag’ (est) database was not considered because it gathers sequences that are both insufficiently accurate and too short (< 800 bp) compared to the length of the HMA genes (> 2200 bp). Among the retrieved sequences, sequences encoding proteins showing orthology to either of AtHMA1, AtHMA2, AtHMA3 or AtHMA4 according to phylogenetic analyses were added to the set of sequences and sequences showing low similarity with our query sequences were discarded. We primarily considered sequence entries corresponding to entire HMA sequences. When an entry only corresponded to a partial HMA sequence, we did not try to establish a full length HMA sequence by assembling that partial sequence with any another partial HMA sequence present in the databases. This was done to avoid assembling overlapping sequence that would actually correspond to different genes. In order to be sure to establish a complete dataset, the identified sequences were translated into protein sequences that were used again as starting query sequences to search for new similar nucleotide sequences in databases, as described above. Only new sequences encoding proteins showing orthology to either of AtHMA1, AtHMA2, AtHMA3 or AtHMA4 were then added to the original set of sequences. This process was repeated until no additional sequence could be added to the set of sequences, considering the sequences available in the databases at the end of April 2011. Some genomic sequences were not considered, either because they were still too fragmented (this was for instance the case for the Lotus japonicus sequences or for the Phoenix dactylifera orthologue to the AtHMA1 sequence) or because they still harboured a great proportion of ambiguous nucleotides (this was, for instance, the situation for the malusXdomesticus sequences).

An important issue was to establish a dataset devoid of duplicates. Sequences were considered to correspond to potential duplicates (as for instance different alleles corresponding to a same locus) on the basis of high sequence identity (> 98%) in the non-coding regions neighbouring the HMA exons. Potential duplicates were removed from the dataset. However, when two or more HMA sequences, showing very high similarity both in coding and non-coding regions were associated in tandem in a single genomic sequence, they were considered as corresponding to distinct genes. This was for instance the case for the three A. halleri HMA4 tandem duplicates, as shown previously [13]. Because of these different filtering steps, sequences reported in Table 1 can be considered as corresponding to distinct paralogues.

A total of 96 HMA sequences were collected (Table 1 and Supplementary Table S1). The dataset comprised sequences from 32 plant species belonging to 15 botanical families (Arecaceae, Asteraceae, Brassicaceae, Caricaceae, Cucurbitaceae, Euphorbiaceae, Fabaceae, Malvaceae, Nicotianeae, Pinaceae, Poaceae, Rosaceae, Salicaceae, Solanaceae and Vitaceae).

One step of the work has been to give names to the HMA protein sequences that were identified. The names started with two letters indicating the plant species from which the sequence originated. Then the three ‘HMA’ letters were used to indicate the gene family. Finally the names ended with either figures (1, 2, 3, 4) or letters (A, B, C, D). We kept their names to the sequences that had already been named. This was the situation for the Arabidopsis, the rice, the N. caerulescens and one of the poplar sequences [1,14,30]. For the previously unnamed sequences, we determined the terminal part of the names after completion of the complete phylogenetic analysis described below. When a clear orthology relationship associated the sequence of a previously unnamed heavy metal ATPase to the sequence of an already named heavy metal ATPase, we put the terminal part of the known sequence at the end of the still unnamed sequence. For instance, the Sorghum bicolor HMA protein sequence extracted from accession number XM_002438908 showed a clear orthology relationship to the OsHMA2 sequence; it was thus named SbHMA2 (Table 1 and Fig. 1). Sometimes two or more previously unnamed sequences showed a clear orthology relationship with one specific already named heavy metal ATPase. For instance, the S. bicolor HMA proteins extracted from accession numbers XM_002459533 and XM_002459534 both showed a clear orthology relationship to the OsHMA3 sequence (Fig. 1). These proteins were then named SbHMA3-1 and SbHMA3-2. When no obvious orthology relationship associated sequences encoding previously unnamed heavy metal ATPases to sequences of already named heavy metal ATPases, the -A, -B, -C or -D letters were appended at the end of the sequence names (for instance, GmHMA-A). As detailed below, this applied to all the HMA sequences corresponding to plant species belonging to other botanical families than Poaceae or Brassicaceae (Fig. 1).

Less than one third of the HMA nucleotide sequences identified as described above corresponded to experimentally validated full-length cDNA clones. The validated sequences were the Arabidopsis and rice ones, as well as a couple of sequences from other plant species such as N. caerulescens, maize, poplar,… (Table 1). Most of the other HMA sequences were collected from on-going genome sequencing projects. They were thus either not annotated or not properly annotated. We thus made the annotation ourselves to obtain the appropriate deduced protein sequences that are given in Supplementary Table S1. In the course of our annotation, we observed that at least four HMA sequences most likely corresponded to pseudogenes. Indeed, these sequences harboured STOP codons within the coding sequence (BrHMA4-2 and GmHMA-D) or displayed only few of the exons (GmHMA-D, PtHMA-B and VvHMA-B). However, in some instances, the genome sequencing was not finished, or the authors who deposited the sequences in the databases only performed partial sequencing. In these situations, we did not consider the incomplete HMA sequences to correspond to pseudogenes and we quoted them as ‘partial’ in Table 1.

From the analysis of the validated set of sequences, we observed that the number of HMA paralogues varied depending on the species considered (Table 1). For instance, most of the species possessed only one HMA1 copy; only Arabidopsis lyrata, barley (Hordeum vulgare) Eutrema parvulum, mustard (Brassica rapa) and soybean (Glycine max) displayed two complete copies. If we now focus on the orthologues of AtHMA2, AtHMA3 or AtHMA4, the situation is more complex. Rice and the strong-spined medick (Medicago truncatula), which genomes are fully sequenced, possess only two paralogues in their genome, while mustard, whose genome is still incompletely sequenced, already displays five paralogues (Table 1). Within the same botanical family, different species display different numbers of HMA copies (compare Z. mays or sorghum to rice for instance, or A. halleri and A. lyrata to A. thaliana in Table 1).

The protein structure of the different HMAs was also highly variable. If we concentrate on the Zn/Cd/Pb/Co P_1B-ATPases other than the HMA1 orthologues for which a full length protein sequence is available, performing a multiple sequence alignment showed that only the first ∼700 amino-acids of the proteins, which harbour the transmembrane domain, display a high level of similarity (data not shown). The cytoplasmic C terminal parts of the proteins appeared to be completely different between orthologues, as well as between paralogues. In particular, they display markedly different lengths between the orthologues, comprising from 58 to 703 amino acids, the average length being ∼290 amino acids (Table 1 and Fig. 2). Intriguingly, some of the Brassicaceae HMAs, the HMA3 ones, display a specific characteristic: their C-terminal part is only ∼60 amino acid long while the C-terminal parts of all the other orthologues to AtHMA2, AtHMA3 or AtHMA4 are much longer (Table 1 and Fig. 2).

A phylogenetic analysis was performed considering the 83 collected HMA proteins for which sequences longer than 550 amino acids were available (Fig. 1). The C-terminal region of the protein could not be considered in this analysis since no sufficient similarity could be detected between the different homologues. Different methods were used including maximum likelihood-, distance- and parsimony-based ones. They all produced nearly identical trees. The very few observed permutations concerned terminal branches of the trees; these permutations did not change the overall organization of the trees (data not shown). This convergence between the results from the different methods warrants the robustness of the phylogenetic analysis.

The phylogenetic analysis revealed that the orthologues of AtHMA1 were all grouped with each other in one clade, whatever their species of origin (Fig. 1). In contrast, the phylogenetic relationship linking the orthologues of AtHMA2, AtHMA3 or AtHMA4 was unusual: the structure of the phylogenetic tree made from these orthologues was strongly dependent from the botanical family to which the sequences could be related to rather than from the different protein isoforms (Fig. 1). Three subgroups were revealed. All the orthologues of AtHMA2, AtHMA3 or AtHMA4 that were identified in Poaceae species constituted a first specific subgroup, within which two separate clades grouping orthologues of OsHMA2 from rice on one side and orthologues of OsHMA3 from rice on the other side could be distinguished. The orthologues of AtHMA2, AtHMA3 or AtHMA4 that were identified in dicotyledonous species were split into two subgroups. One of these subgroups comprised HMAs identified in Brassicaceae species only, and the other one comprised HMAs identified in all the other dicotyledonous species. Within the Brassicaceae subgroup, the HMAs were arranged into three different well-separated clades, one of which contained AtHMA2, the second of which contained AtHMA3 and the third of which contained AtHMA4. In contrast, there was no clear arrangement of the orthologues of AtHMA2, AtHMA3 or AtHMA4 that were identified in non-Brassicaceae dicotyledonous species (Fig. 1); thus these orthologues could not be directly related to either of AtHMA2, AtHMA3, AtHMA4, OsHMA2 or OsHMA3. Finally, three sequences extracted from the monocotyledon P. dactylifera species appear to be scattered in-between the Poaceae sequences and the dicotyledonous sequences.