Introduction
The morphological group of cetrarioid lichens (Parmeliaceae) with erect foliose/subfruticose thallus, marginal apothecia and pycnidia, and production of the Cetraria-type lichenan, contains nearly 150 species in over 20 genera (Randlane et al. Reference Randlane, Saag and Thell1997, Reference Randlane, Saag, Thell and Ahti2010), of which c. 100 species and 17 genera form a monophyletic clade, the so-called ‘cetrarioid core' (Thell et al. Reference Thell, Högnabba, Elix, Feuerer, Kärnefelt, Myllys, Randlane, Saag, Stenroos and Ahti2009; Nelsen et al. Reference Nelsen, Chavez, Sackett-Hermann, Thell, Randlane, Divakar, Rico and Lumbsch2011). Many taxa in the cetrarioid core group are narrowly circumscribed relative to other genera in the Parmeliaceae and so the species and genus delimitations of cetrarioid lichens are still in focus.
The genus Vulpicida J.-E. Mattsson & M. J. Lai, belonging to the cetrarioid core, consists of six species of lichenized fungi that were formerly classified within Cetraria Ach: V. canadensis (Räsänen) J.-E. Mattsson & M. J. Lai, V. juniperinus (L.) J.-E. Mattsson & M. J. Lai, V. pinastri (Scop.) J.-E. Mattsson, V. tilesii (Ach.) J.-E. Mattsson & M. J. Lai, V. tubulosus (Schaer.) J.-E. Mattsson & M. J. Lai and V. viridis (Schwein.) J.-E. Mattsson & M. J. Lai. In Reference Mattsson and Lai1993, the new genus Vulpicida was delimited based on morphology, anatomy, chemistry, distribution and ecology (Mattsson & Lai Reference Mattsson and Lai1993). The genus is distributed in the temperate and arctic regions of the Northern Hemisphere and its species are characterized by a unique set of secondary metabolites, pinastric and vulpinic acids, that are products of the shikimic acid pathway and cause an intense yellow colour of the medulla (Mattsson Reference Mattsson1993). The morphological recognition of the genus is easy due to this impressive character, but its monophyly has not been shown with confidence and the evolutionary relationships between the species have remained unclear.
Vulpicida species have been distinguished using morphological characters and, surprisingly, also distribution. In Northern Europe, three of the six species are known: V. juniperinus, V. pinastri and V. tubulosus (Thell et al. Reference Thell, Ahti, Randlane, Thell and Moberg2011). Vulpicida pinastri is distributed over the Northern Hemisphere whereas the two other species have narrower distribution areas: V. juniperinus occurs in different parts of Eurasia, and V. tubulosus only on alvars of the Baltic Sea islands and in the mountains of Central Europe (Randlane & Saag Reference Randlane and Saag2005). The distribution of the fourth taxon, V. tilesii, is restricted to the arctic and alpine regions of Asia and North America, according to Mattsson (Reference Mattsson1993); however, some herbarium samples from alpine areas in Europe are also labelled as V. tilesii.
Vulpicida juniperinus and V. tubulosus are morphologically very similar, while V. pinastri is easily recognized by being the only sorediate species in the genus. Historically, V. juniperinus was considered to inhabit the branches of Juniperinus communis, and V. tubulosus calciferous ground. Both of them can actually grow on the ground as well as on bark (Mattsson Reference Mattsson1993). Similar overlapping is observed in the morphological characteristics that the authors of the genus used in the species descriptions to distinguish these taxa. In the key for Vulpicida species Mattsson (Reference Mattsson1993: 31–32) used the following characters to distinguish V. juniperinus and V. tubulosus: shape of lobes and lobe margins, colour and distribution. According to him the thallus of V. juniperinus is dorsiventral with thin lobe margins, yellow, and the species occurs in arctic/alpine boreal Eurasia. The lobes of V. tubulosus are terete or dorsiventral with thick margins. Vulpicida tubulosus is greenish yellow and is distributed in the Central European Alps and around the Baltic Sea. In the description of V. juniperinus, Mattsson (Reference Mattsson1993: 38) states that V. juniperinus has only marginal pycnidia on projections, while V. tubulosus can also have laminal and immersed pycnidia. Of anatomical characters, the width of pycnidia was reported as the least overlapping character between V. juniperinus (110 ± 4 µm) and V. tubulosus (90±6 µm), and therefore could be used to separate the two species. In the field, typical V. tubulosus (growing on the ground) is distinguished from typical V. juniperinus (growing on common juniper) by more fruticose growth form with tubular/terete lobes.
Until now, no molecular studies focusing on the phylogeny of Vulpicida species have been published. The previous studies with cetrarioid lichens have shown that the Vulpicida specimens sequenced were related to the taxa from Allocetraria and Cetraria (Mattsson & Wedin Reference Mattsson and Wedin1998; Thell & Miao Reference Thell and Miao1999; Wedin et al. Reference Wedin, Döring and Mattsson1999; Thell et al. Reference Thell, Stenroos, Feuerer, Kärnefelt, Myllys and Hyvönen2002, Reference Thell, Högnabba, Elix, Feuerer, Kärnefelt, Myllys, Randlane, Saag, Stenroos and Ahti2009; Crespo et al. Reference Crespo, Lumbsch, Mattsson, Blanco, Divakar, Articus, Wiklund, Bawingan and Wedin2007; Nelsen et al. Reference Nelsen, Chavez, Sackett-Hermann, Thell, Randlane, Divakar, Rico and Lumbsch2011). However, in these studies, not all species and only few representatives from Vulpicida were included.
In this study we focus on two species: V. juniperinus and V. tubulosus. Distinguishing between these is often very difficult using morphology. Intermediate forms of V. juniperinus and V. tubulosus in Estonia appear to be especially problematic. As a result, we attempt to elucidate the evolutionary relationships of these two taxa using molecular methods.
Material and Methods
Taxon sampling and morphology
The DNA sequences from 79 specimens of eight species were used in this study, 76 of them from five Vulpicida species (Appendix 1). As V. juniperinus and V. tubulosus are the focus of this study, most sequences come from these two taxa (63 specimens together with the morphological intermediates). Cetraria islandica was chosen as the outgroup, based on the recent phylogenetic studies of the cetrarioid core group (Thell et al. Reference Thell, Högnabba, Elix, Feuerer, Kärnefelt, Myllys, Randlane, Saag, Stenroos and Ahti2009; Nelsen et al. Reference Nelsen, Chavez, Sackett-Hermann, Thell, Randlane, Divakar, Rico and Lumbsch2011).
The morphology of V. juniperinus, V. pinastri, V. tilesii and V. tubulosus specimens (74) was observed using a binocular microscope (Olympus SZ51). The main diagnostic characters described by Mattsson (Reference Mattsson1993) were recorded: growth form of thallus, shape of lobes, thickness of lobe margins, presence and type of reproductive structures, and abundance and size of pycnidia. The height and width of an average of five pycnidia per thallus were measured in the case of V. juniperinus and V. tubulosus, but mostly less for V. pinastri and V. tilesii as pycnidia are scarce on these species. To analyze the variance of pycnidial size between clades and morphological types, a nested ANOVA was used in the program Statistica 7.01 (Statsoft Inc.). For this, the height and width of 343 pycnidia were measured on 61 specimens from clades A1, A2 and B. For analysis within morphological types, 29 specimens with typical V. juniperinus or V. tubulosus morphology were selected (15 and 14 specimens accordingly, as shown on the gene tree). Additionally, the colour of thallus and emergence and location of pycnidia were observed on all specimens, but as there was no difference in these characters among the focal taxa of our study, these data are not presented here. From other diagnostic characters, substratum was also recorded.
The content of secondary metabolites was studied in some specimens by means of TLC using solvent system A (Orange et al. Reference Orange, James and White2001); TLC of all specimens was not considered necessary as V. juniperinus and V. tubulosus are chemically identical, according to previous knowledge (Mattsson Reference Mattsson1993; Thell et al. Reference Thell, Stenroos, Feuerer, Kärnefelt, Myllys and Hyvönen2002).
DNA extraction and PCR amplification
The thalli were carefully examined under a stereomicroscope for possible fungal infection. Pieces of the vegetative thallus were used for DNA extraction. The samples were ground in 2 ml microtubes with 2–4 3 mm stainless steel beads using a bead mill (Mixer Mill MM 400, Retsch). Subsequently, total genomic DNA was extracted using High Pure PCR Template Preparation Kit (Roche), according to the manufacturer's instructions with an extra phase separation step using chloroform. The following primers were used: ITS1F (Gardes & Bruns Reference Gardes and Bruns1993) and ITS4 (White et al. Reference White, Bruns, Lee, Taylor, Innis, Gelfand, Sninsky and White1990) for nuITS rDNA (ITS); Mcm7–709for and Mcm7–1348rev (Schmitt et al. Reference Schmitt, Crespo, Divakar, Fankhauser, Herman-Sackett, Kalb, Nelsen, Nelson, Rivas-Plata and Shimp2009) as well as LecMCM7f and LecMCM7r (Leavitt et al. Reference Leavitt, Fankhauser, Leavitt, Porter, Johnson and St. Clair2011) for nuclear Mcm7 gene (Mcm7); mrSSU1 and mrSSU3R (Zoller et al. Reference Zoller, Scheidegger and Sperisen1999) for mitochondrial SSU rDNA (mtSSU). The PCR mix consisted of 12·5 µl PCR Master Mix 2x (Fermentas), containing 0·05 u/µl Taq DNA Polymerase, 4 mM MgCl2, 0·4 mM of each dNTP and reaction buffer, 12·5 pmol of both primers (in 2 µl of water), variable amounts of template and water up to the volume of 25 µl. PCR cycling parameters used for amplifying ITS and mtSSU loci followed Wedin et al. (Reference Wedin, Westberg, Crewe, Tehler and Purvis2009) and for Mcm7 Schmitt et al. (Reference Schmitt, Crespo, Divakar, Fankhauser, Herman-Sackett, Kalb, Nelsen, Nelson, Rivas-Plata and Shimp2009). The PCR products were visualized on 1·5% agarose gels stained with ethidium bromide and purified via SAP/EXO treatment (Fermentas). The same primers were used for sequencing and PCR amplification; both strands of DNA were sequenced using BigDye Terminator v3.1 Ready Reaction Cycle Sequencing Kit (Applied Biosystems). The sequences were run on ABI 3730xl DNA Analyzer (Applied Biosystems). The sequencing procedures were carried out by the DNA Genotyping and Sequencing Core Facility of the Estonian Biocentre and the Institute of Molecular and Cell Biology at the University of Tartu (Estonia). The sequence traces were observed in 4Peaks (mekentosj.com). In some cases, if the sequences of the two strands were not fully complementary, such bases were inspected carefully and usually replaced with ambiguity codes. All sequences were checked using GenBank BLAST.
Sequence alignments and phylogenetic analyses
All regions were aligned separately using the software MAFFT v6.850b (Katoh & Toh Reference Katoh and Toh2008) and MacClade 4.08 (Maddison & Maddison Reference Maddison and Maddison2000). We removed the end of nuSSU, sometimes together with intron, from the beginning of ITS alignment.
The recombination detection program RDP (Martin et al. Reference Martin, Williamson and Posada2005b) was used to scan for possible recombination events in all matrices. The following methods available in this package were employed: RDP (Martin & Rybicki Reference Martin and Rybicki2000), GENECONV (Padidam et al. Reference Padidam, Sawyer and Fauquet1999), Chimaera (Posada & Crandall Reference Posada and Crandall2001), Maxchi (Maynard Smith Reference Maynard Smith1992), Bootscan (Martin et al. Reference Martin, Posada, Crandall and Williamson2005a), SiSscan (Gibbs et al. Reference Gibbs, Armstrong and Gibbs2000), PhylPro (Weiller Reference Weiller1998) and 3Seq (Boni et al. Reference Boni, Posada and Feldman2007).
For phylogenetic inference, we used methods from two paradigms. The single locus gene trees and the concatenated trees were computed following the Bayesian approach (B/MCMC) in MrBayes 3.1.2 (Huelsenbeck et al. Reference Huelsenbeck, Rannala and Masly2000; Ronquist & Huelsenbeck Reference Ronquist and Huelsenbeck2003) and maximum parsimony (MP) with nonparametric bootstrapping in PAUP* 4.0 (Swofford Reference Swofford2002). Bayesian support values can sometimes be overestimates, especially when the tree branches are short. In contrast, bootstrap values can be viewed as lower bounds of support values (Douady et al. Reference Douady, Delsuc, Boucher, Doolittle and Douzery2003). We considered the clades with bootstrap support ≥70% in MP and posterior probabilities ≥95% as strongly supported. The phylogenetic trees were visualized using the program FigTree v1.3.1 (Rambaut Reference Rambaut2009).
Single locus gene trees and concatenated tees
In B/MCMC, independent evolution models were assumed for the matrix partitions. The models were selected using standard AIC (Akaike Information Criterion) in MrModeltest 2.3 (Nylander Reference Nylander2004) and corresponding model form settings were applied in MrBayes. The partitions and selected models were: ITS1 – GTR+G, 5.8S – K80+I, ITS2 – SYM+G, Mcm7 – K80+G, mtSSU – GTR+G.
The gaps were coded as standard characters via modified complex indel coding using the software SeqState (Müller Reference Müller2005, Reference Müller2006). These numeric data were analyzed as separate partitions with the default models for standard characters.
The ITS, Mcm7 and mtSSU regions were first analyzed separately. The numbers of nucleotide positions and the numeric indel codes in the matrices were as follows: 518 and 43 for ITS; 489 and 0 for Mcm7; 839 and 11 for mtSSU. Then the data were concatenated into the three-locus matrix, which consisted of 79 specimens and 1835 nucleotide positions with 54 indel codes. Sequences of all 3 loci were present in the matrix for all 79 specimens.
In all analyses, the MrBayes default priors were used, except for the partition specific rates prior that was set to ‘variable’ (flat Dirichlet). Two simultaneous B/MCMC analyses were run for 10 million (Mcm7, mtSSU) or 12 million generations (ITS, three-locus dataset), both with four chains and starting from random trees. Trees were sampled every 200 generations. The initial 4 million generations were discarded as burn-in from both runs and the majority-rule consensus tree with average branch lengths was calculated for the remaining trees using the sumt command of MrBayes. The average standard deviation of split frequencies between simultaneous runs was 0·005 in three-locus and 0·003 in other analyses; PSRF values all equalled 1·0.
The three-locus dataset was also analyzed by means of MP. Nonparametric bootstrap support (Felsenstein Reference Felsenstein1985) for each clade was estimated based on 1000 replicates, using the heuristic search option with 10 random sequence additions, TBR branch swapping and MulTrees option in effect. The gaps were treated as fifth state.
Results
DNA sequences and recombination
Original sequences were obtained from 79 specimens for ITS, Mcm7 and mtSSU (Appendix 1). Before the ITS1 region, a nuSSU intron was found in all Vulpicida samples except two V. canadensis and two V. tilesii specimens from North America (TIL 03, TIL 04). No credible recombination events were detected, either within or between loci.
Concatenated B/MCMC trees
The B/MCMC consensus tree with posterior probabilities (PP) and maximum parsimony nonparametric bootstrap branch support values (PBS), based on the three-locus matrix (ITS, Mcm7 and mtSSU), is shown in Fig. 1. Starting from the base, three well-supported clades can be distinguished: Allocetraria (PP=100%; PBS=100%), Vulpicida canadensis (PP=100%; PBS=100%) and a clade consisting of four other Vulpicida species, V. juniperinus, V. pinastri, V. tilesii and V. tubulosus (PP=100%; PBS=96%). The latter is divided into two major clades, both consisting of taxa focal to this paper. One of the major clades (PP=100%; PBS=72%) consists of V. pinastri (PP=100%; PBS=95%) and V. juniperinus/ tilesii/ tubulosus groups (clade A; PP=100%; PBS=83%) that in turn includes a few small clades with strong Bayesian and bootstrap support. The second major clade (clade B; PP=100%; PBS=99%) is more distant from V. pinastri and unites the remaining V. juniperinus and V. tubulosus specimens.
Single locus B/MCMC trees
The single locus trees are shown in Figures 2–4. The ITS tree (Fig. 2) was best resolved and supported, followed by Mcm7 (Fig. 3) and the relatively less informative mtSSU (Fig. 4). Vulpicida juniperinus and V. tubulosus were always intermixed. Major topology conflicts above 95% posterior probability level were found between all loci. In the ITS tree, V. pinastri sequences form a strongly supported monophyletic clade (Fig. 2) but, according to Mcm7, some V. pinastri specimens group with a few V. tubulosus and V. tilesii samples (Fig. 3). While in the ITS tree V. juniperinus and V. tubulosus form two distinct intermixed clades, in the Mcm7 analysis most sequences of these taxa are found in one clade with a smaller (also mixed) subclade that has no match in the ITS tree. In the Mcm7 tree, Allocetraria appears closer to the V. juniperinus/ pinastri/ tilesii/ tubulosus group than V. canadensis, while in the ITS tree V. canadensis groups with one of the V. juniperinus/ tubulosus clades. According to mtSSU, all Vulpicida specimens are united in one clade with low support but the inner structure of that is different from the ITS and Mcm7 trees.
Morphological characters
By assigning species epithets to the specimens, we tried to reflect the morphological traits of the thalli as much as possible and therefore intermediate specimens were marked with more than one epithet in the order of their affinity (Appendix 1).
The states of the recorded morphological characters are given in Fig. 1. Both V. juniperinus and V. tubulosus were found growing on bark as well as on the ground. The V. juniperinus growth form is foliose to subfruticose, and V. tubulosus subfruticose to fruticose. Few thalli had terete lobes only; most thalli of V. tubulosus had a mix of terete/ radial and flat/ dorsiventral lobes (Fig. 5.). Apothecia were found on only a few thalli. Vulpicida. pinastri is distinguished by having soralia. Pycnidia were absent only on one studied thallus. No differences in thallus chemistry between the tested specimens were found by means of TLC.
The height and width measurements of pycnidia were as follows [(min) average±standard deviation (max)]. Height in clades A1+A2: (48) 118±26 (172) µm; clade B: (72) 117±24 (168) µm; typical V. juniperinus: (72) 125±24 (168) µm; typical V. tubulosus: (80) 116±24 (160) µm. Width in clades A1+A2: (80) 112±19 (160) µm; clade B: (80) 105±20 (164) µm; typical V. juniperinus: (80) 115±18 (164) µm; typical V. tubulosus: (72) 102±19 (144) µm.
Nested ANOVA analysis was carried out to compare the sizes of pycnidia between: 1) clades A1+A2 and clade B, and 2) morphologically typical specimens of V. juniperinus and V. tubulosus that were first determined without measurements of pycnidia. There was significant variation of the height (F 59282= 1·999, P<0·001) and width (F 59282= 1·440, P=0·027) of pycnidia among specimens, but only the variance of pycnidial width was significant among clades A and B (F 1282=10·810, P<0·001). There was also significant variation in the height (F 1144=4·893, P<0·029) and width (F 1144=18·460, P<0·001) of pycnidia among morphological groups.
Discussion
Despite the conflicts between loci, our concatenated three-locus tree is surprisingly well resolved and supported, clearly better than any single locus tree. However, we cannot be sure if this method handles the conflicts correctly. The results of the analysis are dominated by the strong signal from ITS. It shows that V. juniperinus and V. tubulosus specimens are divided into two clearly distinguished groups by the DNA characters. One of these groups (clade A) is more closely related to the monophyletic V. pinastri than to the remaining V. juniperinus and V. tubulosus sequences (clade B). The two morphospecies are mixed between the clades and V. juniperinus and V. tubulosus appear polyphyletic in the concatenated, as well as in single locus, analyses. A tendency can be seen that more specimens with the typical V. juniperinus morphology are located in clade A than in clade B, which includes more morphologically characteristic V. tubulosus specimens. However, representatives of both morphospecies are clearly present in both clades. Moreover, a few morphological V. juniperinus and V. tubulosus specimens have identical sequences in all three loci. Thus, the main outcome of our study is that currently accepted V. juniperinus and V. tubulosus cannot be separated using the DNA characters of our studied loci.
We found no distinctive morphological characters for clades A and B except the width of pycnidia. However, this difference is small and the large variability and overlap of the measurements between the groups make this character impractical. Also, in many cases the width of pycnidia is clearly incompatible with the rest of the characters, in respect of the morphological species concepts currently applied. As the width and height of pycnidia are significantly different between morphospecies, the significant difference in the width of pycnidia between the two clades could be partly due to uneven dispersal of morphological types between clades.
There seems to be no significant correlation between the clades of the concatenated trees and geography concerning V. juniperinus and V. tubulosus specimens from Europe. Vulpicida juniperinus specimens from Japan form a strongly supported monophyletic subgroup. Between the subclades of V. pinastri, no geographic distinction was detected.
Sequences of V. pinastri form a strongly supported clade. This group is easily separable from others also by morphology (soralia on thalli). The few (four) V. tilesii specimens sampled do not form a monophyletic clade in the concatenated or single locus trees but cluster together with V. juniperinus and V. tubulosus specimens in clade A. The position of V. tilesii needs to be clarified in the future using a greater number of specimens.
Morphologically similar but genetically distinct species are sometimes hidden under a single name. Such phylogenetic species are frequently referred to as ‘cryptic’ (Funk & Omland Reference Funk and Omland2003; Crespo & Pérez-Ortega Reference Crespo and Pérez-Ortega2009). The discordance between traditional morphological species concepts and DNA based phylogeny is not uncommon (Lumbsch & Leavitt Reference Lumbsch and Leavitt2011). One could conclude, based on our ITS and concatenated trees, that V. pinastri splits two morphologically indistinguishable cryptic species that correspond to clades V. juniperinus/ tubulosus A and V. juniperinus/ tubulosus B. However, this might not necessarily be the case. Relatively short branches in the phylogenetic trees and little differences in other characters indicate that we might be dealing with a young species complex consisting of V. juniperinus and V. tubulosus together with V. tilesii and V. pinastri. In young diverging taxa, a common cause of the apparent polyphyly in gene trees is the incomplete lineage sorting (Grube & Kroken Reference Grube and Kroken2000; Taylor et al. Reference Taylor, Jacobson, Kroken, Kasuga, Geiser, Hibbett and Fisher2000; Funk & Omland Reference Funk and Omland2003; Knowles & Carstens Reference Knowles and Carstens2007). Also, when there is insufficient time for the fixation of gene lineages, the random division of allele copies during speciation can lead to different gene histories between loci (Hudson & Coyne Reference Hudson and Coyne2002), and thus to the incongruence between the single locus gene trees that we found in this study.
One of the major biological reasons for polyphyly of a species is considered to be recombination, which can be detected using numerous algorithms. We found no credible recombination events within or between loci; however, we cannot completely rule out their presence in our matrix, as different recombination detection methods show distinct performances depending on the amount of recombination, genetic diversity and other properties of the sequence data (Posada & Crandall Reference Posada and Crandall2001).
As the division of V. juniperinus and V. tubulosus sequences into the clades A and B is only present in the ITS tree, not Mcm7 and mtSSU one-locus trees, one could consider this to be the result of paralogy: two slightly different copies of the same locus in one genome. Even though we cannot exclude this option entirely without applying cloning, we are not seeing clear signs of paralogy in our study. From some V. juniperinus and V. tubulosus specimens, the ITS region was independently amplified and sequenced twice with different template dilutions (some also extracted twice) resulting in identical sequences (4 from both clades A and B). From few other repeated PCRs the sequences of parasitic fungi in addition to the ones of Vulpicida were obtained, as identified by BLAST, but never any possible paralogs. The variable sites in the ITS matrix distinguishing clades A and B are all based on high quality sequence traces and none of those show traces of the alternate version of the locus.
Besides paralogy or recombination, it is also possible that the clades A and B really are in an early stage of speciation that is reflected only in ITS and not yet in the other two markers because of different molecular evolution rates between loci. But in this case we still would not call them species before other loci support it. We share the view of Grube & Kroken (Reference Grube and Kroken2000) that the species should be separated based on single locus only if the clades clearly correlate with phenotypic characters or biogeographic distribution.
The results of our analyses do not support the two traditional species separated by morphology. As the samples of these species are intermixed in all gene trees, distinguishing V. juniperinus and V. tubulosus according to the current descriptions is not justified. We do not propose nomenclatural changes in this paper because the phylogenetic situation in the genus Vulpicida is complex, and to overcome the conflicts between the one-locus trees, more loci and more phylogenetic methods should be employed. Before we can delimit any taxa with confidence, other Vulpicida species, especially V. tilesii, need to be sufficiently sampled.
We thank the collectors of the specimens, especially Emma Sandler Berlin and Ede Leppik, and the curators of GZU, LD, MAF and UPS for the loans. The study was financially supported by the Estonian Science Foundation (grants JD173, 7470 and 9109), and the European Union through the European Regional Development Fund (Center of Excellence FIBIR). Some of the phylogenetic analyses were carried out in the High Performance Computing Center of the University of Tartu.