Proteomic similarity of the Littorinid snails in the evolutionary context

Ethics statement

Periwinkles were collected from wild populations (Table 1). Periwinkles are not endangered or protected species in the study regions, so no special permission for their collection was required. Administrations of the Swire Institute of Marine Science, Hong Kong University (Hong Kong, China), the Arctic University of Norway University of Tromsø (Norway), the Educational and research station “Belomorskaia” of St. Petersburg State University (Russia), the Scotland Association of Marine Science (Oban, Scotland, UK) and Interuniversity Institute for Marine Sciences (Eilat, Israel) were informed about the snails sampling. Collection and animal numbers were approved by these authorities. No special approval could be obtained for sample collection at the site near Cancale (France) because it is not an area of any national park or a private territory, so there was no appropriate authority to apply for permission.

Sample collection and preparation

Adult individuals of littorinid snails Littorina littorea (Linnaeus, 1758), Littorina fabalis (Turton, 1825), Littorina obtusata (Linnaeus, 1758), Littorina arcana Hannaford Ellis, 1978, Littorina compressa Jeffreys, 1865, Littorina saxatilis (Olivi, 1792), Littoraria ardouiniana (Heude, 1885), Littoraria “melanostoma E Asia” (is unnamed, see (Reid, Dyal & Williams, 2010)), Echinolittorina millegrana (Philippi, 1848) and Echinolittorina marisrubri Reid, 2007, were collected from wild populations in six geographic locations (Table 1).

After collection animals were transported to the laboratory. The Littorina snails from the Northern Atlantic were kept in aerated moist containers at 8 °C and rinsed with salt water once a day. Snails were acclimated to these standard conditions for no less than 2 days and no longer than 1 week before dissection. The Echinolittorina and Littoraria snails from the Red and East-China Sea coasts, respectively, were kept in dry containers at room temperature and rinsed with salt water once a day.

For sample preparation snails were dissected under MBS-10 binocular microscope at 79-849 magnifications to identify species, sex and possible trematode or other heavy parasitic infection and to separate body parts (head, foot and penis). The mollusks identification was conducted according to the shell form, sculpture and pigmentation, as well as male and female reproductive system anatomy (Reid, 1986, 1996, 2007; Granovitch et al., 2008). The correctness of the tropical species identification was confirmed by Dr. David G. Reid personally. The identification of the species within the Littorina “saxatilis”- and “obtusata”-groups (the most challenging part) was performed as described previously (Reid, 1996; Granovitch et al., 2008). More details can be found in the File S1. Only the reliably identified mature individuals with well-developed reproductive system and free of trematode or other obvious infection were used in further analysis.

Tissues of up to 20 animals (collected from different parts of intertidal area within the same location) were pooled separately for different body parts (foot, head, penis) and sexes; and homogenized in lysis buffer (7 M urea, 2 M thiourea, 4% CHAPS, 25 mm Tris, pH 8.2) using Mixer Mill MM 400 (Retsch). Particles were sedimented by centrifugation 12,000g, 15 min, 4 °C and supernatants were frozen at −80 °C until use.

Proteomic analysis strategy

The “biological noise of the system” due to individual variations and microniche adaptation was leveled by pooling, the short acclimation procedure and including biological replicates into the analysis (rev. in Karp et al. (2005)). Samples of somatic tissues of males and females were used as biological replicates, as a previous study showed no significant sex effect at the qualitative level (Maltseva et al., 2016). Head and foot tissues used for analysis in this study are less sensitive to transient environmental impacts than for example, gill as was demonstrated earlier in a study on salinity stress (Muraeva et al., 2016). Technical replicates (at least two for every sample) were averaged to reduce “technical noise” (Karp et al., 2005). Altogether this strategy ensures the analysis of the most stable and reproducible proteome, reducing accidental impacts.

In general, the qualitative strategy of a proteomic study aims for the characterization of the set of proteins present in a sample. This was successfully applied to characterize proteomes of different compartments (organelles, tissues, organs) and certain classes of proteins etc., being an integrative part of structural biology. The quantitative approach estimating proteins abundances is a part of functional biology, and is more informative for studies of physiological changes, adaptive mechanisms, etc. (rev. in Messana et al. (2013)). We chose the qualitative approach to compare basic proteomes between species since it is less amenable to transient physiological and environmental impacts. Our approach principally resembles the blind whole-genome assaying methods like AFLP, RFLP, RAPD, etc. in the sense of that the identity of analyzed proteomic markers (protein spots) is unknown. The main difference is taking into account the “working part of the living molecular machinery”, that is, the proteome, while abovementioned methods engage DNA-regions, both coding and noncoding. We separated proteins in two-dimensional gel electrophoresis, and ultimately, estimated the number of common and unique (new or altered due to amino acid substitutions, modifications, splicing patterns, etc.) protein signals in single species or species groups in comparison (more details on our approach and rationale of proteomic vs. genetic data comparison can be found in the File S1).

We followed a gel-based methodology because it has, among others, an advantage of whole-protein analysis and consequently saves combinations of different post-translational modifications as individual signals (Arentz et al., 2015). We used differential 2D-electrophoresis (DIGE (Ünlü, Morgan & Minden, 1997)) because it allows the simultaneous analysis of up to three samples within the same gel (allowing the reliable detection of even minor structural changes) and provides high sensitivity of protein detection (down to 150–500 pg of a single protein) (Lilley & Friedman, 2004). Qualitative gel analysis and spot matching were carried out using PDQuest Advanced 8.0.1 software (BioRad, Hercules, CA, USA) by making whole experiment master gel, allowing for the warping of individual gel scans (Berth et al., 2007). Spots were considered reliably detected if they were detected in at least two technical replicates of the same sample, or in one technical replicate in at least two different samples. We have taken into account potential post-translational modifications and/or splice variants by identifying independent signals, because every particular modification/splice form possesses its own features and functions, and appearance of a new form/modification represents an important evolutionary event, comparable to amino acid(s) substitution or even the emergence of a new protein. In total, 796 spots were used for inter-species comparison of proteomes. Protein spots differing between species and available for excision after Coomassie gel staining were subjected to tandem mass spectrometry MS/MS-analysis for identification following standard “bottom up” protocol. Basically, in this method proteins are digested into peptides prior to MS/MS-analysis, molecular masses of resulting peptides and their fragments (after fragmentation in the collision cell of mass-spectrometer during MS/MS-analysis) are determined and compared with MS-spectra predicted from a database by software (see File S1 for more details). The success or failure in protein identification did not influence the results of inter-specific comparisons, and aided only in their interpretation (more details on experimental procedures can be found in File S1).

Statistical analysis

All analyses were run using R (R Core Team, 2015).

Phylogenetic analysis

We used sequence data from the previous extensive study on littorinid phylogeny (Reid, Dyal & Williams, 2012) to produce a tree only for those species that are included in the present study. Fragments of 28S rRNA, 12S rRNA and cytochrome oxidase subunit I (COI) were obtained from NCBI; information about sequences used is available in File S1. Sequences were concatenated with different genes allowed to evolve independently. Phylogenetic analysis was performed using Bayesian inference, MrBayes 3.2.6 (Ronquist & Huelsenbeck, 2003) with the same parameters as in the original study (Reid, Dyal & Williams, 2012): fragments of the three genes were used as a concatenated sequence, where individual genes were unlinked to evolve independently; nucleotide substitution model was GTR + G + I; analysis was performed as two independent runs, five chains in each (four heated and one cold; the first 25% samples from the cold chain were discarded) for 25,000,000 generations with a sample frequency of 1,000, print frequency of 1,000 and diagnostics calculated every 1,000 generations. The convergence between two runs was tested by comparison of statistical parameters in the Tracer Software (http://tree.bio.ed.ac.uk/software/tracer/).

Comparison of phylogenetic and proteomic trees

Prior to comparison, phylogenetic and NJ-based proteomic trees were made ultrametric using non-negative least squares with the phangorn package (Schliep, 2011). The Robinson–Foulds distance (both raw and normalized) between unrooted trees was computed to assess the total number of partitions which are implied by one of the trees, but not the other (Robinson & Foulds, 1981). Compositional similarity of proteomic and phylogenetic trees was evaluated using the Fowlkes–Mallows index (Fowlkes & Mallows, 1983). Correspondence between inter-species distances on the trees was measured using cophenetic correlation (Sokal & Rohlf, 1962). In addition, cophenetic distances on the proteomic and phylogenetic trees were expressed as fractions of maximal tree distance to allow direct comparisons (the matrices of transformed distances are presented in Fig. 1).

General trees’ topologies

We applied two algorithms for proteomes clustering (NJ and UPGMA) which both gave similar results (Figs. 1 and 2; File S2). We deduced the overall similarity of species proteomes from their clustering based on presence/absence of peptides in consensus proteomes (Fig. 1A; File S2). These results were compared with a phylogenetic tree obtained by Bayesian inference based on three molecular markers used earlier by Reid, Dyal & Williams (2012) (Fig. 1B; File S2). The general trees’ topologies were similar: the first split was between genera—Littoraria, Echinolittorina and Littorina (however, the genera were split in different order at the proteomic and genetic-based trees); then within Littorina subgenera (L. (Littorina) littorea and L. Neritrema spp. separated), and finally the “saxatilis”-group split from the “obtusata”-group. This was confirmed by the low Robinson–Foulds distance (raw RF = 2, normalized RF = 0.43), high cophenetic correlation coefficient (0.80), and by the distribution of the Fowlkes–Mallows index, which indicated significant compositional similarity for all numbers of clusters, except nine (Fig. 1C). The discrepancy between the trees with nine clusters came from the “saxatilis”—group species cluster, where three species grouped differently (Fig. 1): L. compressa was a sister-group to the pair L. saxatilis/L. arcana in the proteomic tree and L. saxatilis was a sister taxon to the pair L. compressa/L. arcana in the genetic one. Noteworthy, clustering in the proteomic tree was supported by high bootstrap values, unlike low support (0.66) of the “saxatilis”—cluster in the DNA-tree.

Interspecies cophenetic distances

On the whole, proteomic distances in general exceeded genetic distances (the full matrices of distances, both genetic and proteomic, are available in Fig. 1D). The most robust differences were as follows. Relationships in the “saxatilis” group species cluster: L. arcana and L. compressa were considered phylogenetically close, but rather distant at the proteomic level (the proteomic and genetic distances between L. arcana and L. compressa were 0.133 and 0.017, respectively); on the contrary, L. arcana was proteomically very similar to L. saxatilis in contrast to the genetic distance (the proteomic and genetic distances between L. arcana and L. saxatilis were 0.021 and 0.034, respectively); the differences were even more profound if UPGMA clusterization algorithm was applied (see File S2). Distances in the “obtusata” group species cluster: L. fabalis and L. obtusata were identified as phylogenetically very close but strongly diverged at the proteomic level (the proteomic and genetic distances between L. fabalis and L. obtusata were 0.285 and 0.013, respectively). Their proteomic cophenetic distance is comparable with one between genetically distant Littoraria ardouiniana and L. melanostoma (Fig. 1A, Littoraria species formed common cluster just one step earlier than L. obtusata/L. fabalis, Table 2; Table S2). The Littoraria genus case. Similarly, proteomic distances were larger than genetic ones for two Echinolittorina species and two Littorina subgenera (Littorina and Neritrema). At the same time, two Littoraria species showed an unusual outcome: being phylogenetically quite distant, they appeared a bit closer to each other at the proteomic tree (the proteomic and genetic distances between L. ardouiniana and L. melanostoma were 0.313 and 0.324, respectively), which might be interpreted as a sign of physiological and molecular parallelism.

Effect of geography, sex and body part

To assess the effect of geography, sex and body part on the similarity of proteome composition, we calculated the extended proteomic tree of samples (Fig. 2; File S2). The general topology of this tree was similar to the consensus one. Almost all species formed clear species-clusters, which supports the potential taxonomic use of proteomics. Nevertheless, there was an exception: the pair L. saxatilis/L. arcana formed a mixed cluster, where the factor “body part” (and sometimes “geographic location”) affected proteome stronger than “species”. This cannot be due to male species-misidentification, as somatic samples of those males grouped together with unmistakably identified females of these two species. The first within-cluster division splits the UK samples of both species from the mainland ones, suggesting that the geographic variation surpassed the interspecific variability, and a parallel pattern is observed in the proteomes of these two species. Remarkably, distinct UK-clades were also present within the species-clusters of L. obtusata, L. fabalis and L. compressa.

Cluster stability of samples and cluster markers

To further assess the proteome species-specificity we analyzed the stability of sample clustering (Hennig, 2007). Splitting the bootstrapped tree of samples into 10 clusters (the number of species analyzed) resulted mainly in same-species clusters with high recovery probability (Table 2; Table S2), with the exception of the L. arcana/L. saxatilis pair which formed two multi-species clusters, corresponding to European mainland and UK locations (both with maximal recover probability). All species branches were supported by the presence of highly reliable (i.e., maximum specificity and sensitivity) protein markers, except for L. saxatilis (Table 2), which appeared poorly distinguishable from L. arcana at the proteomic level. For L. arcana there were only two specific markers. Moreover, these two markers represent two forms of the same enzyme: the peptidyl-prolyl cis-trans isomerase, which additionally reflects the physiological similarity of these two species.

Markers analysis

We succeeded in the MS/MS-identification of several species (cluster) markers (of 796 protein spots analyzed 128 were successfully annotated, Tables S1 and S3). Expectedly, the most efficient identification was based on a L. saxatilis database (LSD) search (Canbäck et al., 2012) and gave positive results mostly in the case of L. saxatilis and L. arcana (~50% efficacy) in comparison to other species (~4–20% efficacy, Table S3). A variety of house-keeping proteins were usually revealed as markers of cryptic groups (“saxatilis” and “obtusata”) and genera (Littoraria and Echinolittorina species pairs), such as protein disulfide isomerase, calreticulin, tropomyosin, troponin T, radular globin, etc. Metabolic enzymes and stress related proteins (both potentially adaptively important), like arginine kinase, fructose-bisphosphate aldolase, malate dehydrogenase, peptidyl-prolyl cis-trans isomerase (PPIase), glutathione-S-transferase, etc., were often identified as species-branch supporting markers (the latter was also the case for L. saxatilis/L. arcana pair vs. L. compressa). The unusual pattern was found for the Littoraria species pair, where species with quite long independent evolutionary history matched by core metabolic enzymes (unlike, e.g., “saxatilis” group), but differed by globins (which were common within either the “saxatilis” or “obtusata” group and within the Echinolittorina species pair).

Similar proteomes in deeply diverged species: a possible case of the whole-proteome convergence/parallelism

In general, the proteomic distances in the trees exceeded the genetic ones. This was rather expectable based on the inclusion of fast evolving proteins; for example, arginine kinases (present as several forms in periwinkles proteomes (Maltseva et al., 2016)) were among annotated proteins, and these enzymes are well known for their non-conservative nature (Uda et al., 2006). However, there was one remarkable exception from the “large proteomic distances”-rule. Two deeply diverged Littoraria species appeared in the proteomic tree closer to each other than in the genetic one. This result may be interpreted as a molecular parallelism. Although there are many known examples of molecular convergence/parallelism (reviewed in Storz (2016)), usually this phenomenon is analyzed at the level of individual proteins of particular families, not the whole proteome. There are strong theoretical arguments for molecular homoplasy resulting from neutral and not adaptation-driven phenomena (Zou & Zhang, 2015), for example, a set of structural constraints, limiting protein changes to certain trajectories (reviewed in Storz (2016); Starr & Thornton (2016)). Interestingly, in the present study another tropical species pair E. marisrubri/E. millegrana (their estimated divergence-time is slightly shorter than that of the Littoraria species, ~35 Ma and ~48 Ma, respectively) (Reid, Dyal & Williams, 2012) showed no sign of proteomic parallelism, though both species pairs live under extreme conditions of thermal stress, hypoxia and desiccation. Nevertheless, adaptation to a rather special biotope in mangroves, sympatrically inhabited by L. ardouiniana/L. melanostoma, could be another possible explanation, as both species are ecologically very similar in terms of host-tree species distribution, diet and activity (Lee, Williams & Hyde, 2001; Lee & Williams, 2002a). Altogether, the example of the two Littoraria species shows that similar proteomes can emerge from divergent genomes.

Interestingly, globin proteins deviate from the observed pattern of whole-proteomic similarity, showing differences between the two Littoraria species. This may appear quite puzzling; for comparison, no globin-related difference was observed between the strongly diverged Echinolittorina species. A possible explanation for the globin divergence in Littoraria species may be differences in their vertical distribution and expected different aerial-aquatic respiration VO₂-ratio (McMahon, 1988; Lee & Williams, 2002b). Similarly, the “obtusata”- and “saxatilis”-groups differ in their preferences to littoral zone level (Granovitch et al., 2008, 2004). Accordingly, they demonstrate different VO₂-ratio (McMahon, 1988) and different globins, though these species are much more phylogenetically close compared to the Echinolittorina species pair.

Distant proteomes in phylogenetically close species with different habitat preferences

The species pair of the “obtusata” group—Littorina obtusata and L. fabalis shows a pattern that is opposite to the one described above for Littoraria. In the Littorina “obtusata” group we observed strongly diverged proteomes, apparently originating from very akin genomes in morphologically similar species. These two species are phylogenetically very close (closer than any other in the Littorinidae family; estimated divergence time is less 1 Ma) (Reid, Dyal & Williams, 2012) and were shown to share some mitochondrial haplotypes (Kemppainen et al., 2009) (and there are signs of limited gene flow between these two species in Portugal (Costa et al., in press)). Here, the comparative proteomic analysis revealed robust differences between L. obtusata and L. fabalis. The distance between them highly exceeds those within “saxatilis”-group. The discrepancy between genetic and proteomic similarities can be an indirect evidence of ecological speciation. L. obtusata and L. fabalis demonstrate different preferences in their microhabitat distribution: L. fabalis occupies lower intertidal and subtidal levels, while L. obtusata prefers the middle littoral zone; at those levels, they are associated with different fucoids (Reid, 1996; Granovitch et al., 2008, 2004). One can hypothesize that partitioning of ecological niches drove divergence at the physiological level, finally reflected in very different proteome patterns in these two species. The observed striking proteomic differences were achieved by these two species in relatively short evolutionary time, which implies a high evolvability of the system. However, the exact nature of proteomic differences is unknown (whether it is related to genes duplication, amino acid substitution, gene expression pattern shift, alternative splicing or some posttranslational modification), and at least part of the differences can be due to phenotypic plasticity. Thus, these species present an interesting system to further investigate mechanisms of rapid proteome divergence.

Distant proteomes in phylogenetically close species with similar habitat preferences

The third example, from the Littorina “saxatilis” group, suggests that different proteomes can be generated by similar genomes even if corresponding species populate very similar biotope. The three species in this group (L. arcana, L. compressa and L. saxatilis) have a high degree of whole genome similarity, as was demonstrated by comparative genomic hybridization (Panova et al., 2014). Earlier molecular phylogenies were inconclusive regarding these species relationships (Canbäck et al., 2012; Reid, 1996; Reid, Dyal & Williams, 2012) and different molecular markers clustered the three species of the “saxatilis” group in all possible combinations (Panova et al., 2014; Reid, Dyal & Williams, 2012; Wilding, Grahame & Mill, 2000a, 2000b; Small & Gosling, 2000; Knight & Ward, 1991), which is one line of evidence for their close phylogenetic proximity and a recent evolutionary divergence (divergence time is estimated as ~1–1.5 Ma) (Reid, Dyal & Williams, 2012).

In our proteomic tree, the first divergence event splits L. compressa from a tight alliance of L. saxatilis/L. arcana (Fig. 1A). We found 11 species-specific markers for L. compressa and a similar number of markers specific for the pair L. saxatilis/L. arcana, but only two markers specific for L. arcana and none for L. saxatilis (Table 2). These results are congruent with the fact that both L. saxatilis and L. arcana demonstrate parallel proteomic shifts along the vertical shore gradient while L. compressa does not (Maltseva et al., 2016). The separation of L. compressa at the proteomic level partially resembles the foregoing case of L. fabalis and L. obtusata in a sense that it occurs between very closely related species. However, unlike the latter pair, the proteome divergence is not easily explained by ecology in this case. Both L. arcana and L. compressa are very similar to L. saxatilis: their ranges and habitat distribution overlap (Reid, 1996; Granovitch et al., 2013). Although L. compressa was not recorded from fucoid-free substrates, unlike two other species, it invariably lives in sympatry with L. saxatilis. All three species are often found together on the same fucoids and rocks within an intertidal area, and there are so far no data suggesting stronger differences in ecological preferences between L. compressa and L. saxatilis than between L. saxatilis and L. arcana. It is still possible that L. compressa occupies very specific microhabitats to be elucidated in future studies. Alternatively, some hypothetical switch could have occurred in the molecular machinery functioning of L. compressa. The last explanation suggests that speciation events could be related not only to environmental shifts or niche expansion but also to changes in an organism’s patterns of interaction with the same (or very similar) environment. Among protein markers specific to L. compressa we identified arginine kinase, fructose-bisphosphate aldolase and glutathione-s-transferase. These enzymes are potentially related to stress adaptation (e.g., isolation-related hypoxia during low tide; in details discussed in Maltseva et al. (2016)); and such changes are compatible with any of the two proposed explanations.

Strinkingly similar proteomes in phylogenetically close species

In contrast to the L. compressa case discussed above, the pair L. saxatilis/L. arcana showed high proteome similarity. In general, this similarity is expected pattern given the short divergence time and overlapping habitats. What is surprising is that both “geographic location” and “body part” showed stronger impacts on the proteome than “species”; that is, the first split occurs between UK and other regions (further discussed below), then between “foot + head” and “penis” samples, and then between the species (Fig. 2). The misidentification of morphologically similar males of L. saxatilis and L. arcana could not explain this pattern, because female “foot + head” samples cluster together with male “foot + head” samples; and females of the two species are unambiguously distinguishable.

L. saxatilis and L. arcana have strikingly similar proteomes (at least in adults and the analyzed body parts), even though they have different reproductive strategies (L. arcana is oviparous while L. saxatilis is ovoviviparous) and demonstrate some differences in distribution, morphology and genetics (Panova et al., 2014; Reid, Dyal & Williams, 2012; Wilding, Grahame & Mill, 2000a; Granovitch et al., 2008; Knight & Ward, 1991; Mikhailova et al., 2009). Given their genomic, morphological, ecological and now even proteomic similarity, this pair of species raises intriguing evolutionary questions. (1) What mechanisms maintain this similarity across their range? (2) What kind of barrier separates them from each other given indirect evidence of gene flow between them (Granovitch et al., 2013; Mikhailova et al., 2009)? (3) What evolutionary forces drove their divergence in the first place?

Only two species-specific protein markers were detected in our analysis for L. arcana (and none for L. saxatilis). Both L. arcana-specific markers (and one marker of L. compressa) proved to belong to PPIase class of enzymes, the cyclophilins (Cyp) family. These ubiquitously distributed enzymes are known to facilitate folding of certain proteins and thus regulate diverse cellular processes (rev. in Galat & Metcalfe (1995); Fischer & Aumüller (2003); Lu et al. (2007)). In mollusks, Cyps are expressed in hemocytes and gonads and involved in immune response (Song et al., 2009; Ji et al., 2013). They also participate in the processes of biomineralization during shell and radula formation (Jackson et al., 2009; Marie et al., 2013; Nemoto et al., 2012). These enzymes are believed to be secreted to spatially organize matrix proteins, and thus to control crystallization. Hypothetically, structural differences in Cyps between closely related species could cause tiny changes in this process resulting in interspecies radula variation, which was documented for L. arcana and L. compressa (Reid, 1996). In our analyses we used head tissues including odontophore and radular invagination, where radula synthesis occurs. Moreover, Cyp of L. saxatilis is expected to be catalytically active as there are no substitutions in catalytically important amino acids (based on Cyp-transcript sequence present in LSD transcriptome (Canbäck et al., 2012)). Further studies are needed to examine the precise compartment of expression and functioning of the identified Cyps.

Legacy of post-glacial history in proteome divergence

Overall, samples from different localities grouped by species implying a tight link between genome and proteome and a minor influence of physical factors, locality or season. Nevertheless, the same clustering pattern repetitively appeared within every Neritrema species: UK samples formed a distinct clade (Fig. 2), while continental samples clustered more often according to the tissue type. In the pair L. arcana/L. saxatilis separation between UK and continental samples even precedes the separation of the species (Fig. 2; Table S2). L. saxatilis UK populations have been previously shown genetically distinct from the mainland using mitochondrial markers (Panova et al., 2011; Doellman et al., 2011) and it has been suggested that UK L. saxatilis persisted in the UK within some local refugium during the last glacial maximum with following postglacial expansion on the British Islands, but not to the mainland. European mainland populations of L. saxatilis were presumably recolonized from other refugia than the UK (Panova et al., 2011; Doellman et al., 2011). The existence of a UK refugium has been also suggested for L. compressa and L. arcana (Doellman et al., 2011). Our findings are well in line with such scenario. Moreover, both “obtusata” species have possibly also survived in a hypothetical UK refugium, as a separate UK-clade was also detected in these two species. Recently, a new phylogeographic study suggested the UK as a possible glacial refugium for L. fabalis and L. obtusata (Sotelo et al., in press). Overall, we observed a pattern of geographic isolation at the proteome level, which is very similar to effect of the recent post-glacial history observed in mtDNA. It would be interesting to test this effect further by proteomic analyses of other populations known to have survived in separate glacial refugia, such as L. saxatilis, L. obtusata and L. fabalis in the Iberian peninsula, or L. saxatilis and L. obtusata in North America.

Qualitative proteomics is not reliable for intergenera analysis

When qualitative data of different genera are put into comparative analysis, proteomics cannot resolve their relationship (unlike the interspecies analysis). This occurs because too many proteins are already changed and too few are conserved, yielding a low number of common “anchoring” proteins between strongly diverged genera and poor grouping of the genera with comparably high binary dissimilarity. This is well illustrated by short genera branches relative to species ones (Figs. 1 and 2). Thus, although in our study genera clustered differently than in DNA-markers-based tree, this result is not sufficient to claim a conflict between the two approaches in this case.