Molecular Phylogeny and Phylogeography of Potentilla multifida L. agg. (Rosaceae) in Northern Eurasia with Special Focus on Two Rare and Critically Endangered Endemic Species, P. volgarica and P. eversmanniana

The results of a molecular genetic study of Potentilla multifida agg. using two plastid markers (ndhC-trnV and psbA-trnH) and a nuclear ITS marker suggested that this group comprises a number of relatively young and incompletely differentiated species widely distributed in Northern Eurasia. The sequences were analyzed using tree-based (maximum likelihood) and network-based (statistical parsimony network) approaches. The plastid data suggested incomplete lineage sorting, characteristic of the group as a whole. The nuclear ITS results demonstrated quite a different pattern, with mostly conspecific accessions shaping monophyletic clades. The majority of the Potentilla sect. Multifidae species studied possess few, usually closely related plastid haplotypes, or are even monomorphic. In contrast, P. volgarica, a narrow endemic from the Volga River valley, presents plastid haplotypes belonging to two distantly related groups. Such a pattern of genetic diversity in P. volgarica may be explained by a long persistence of the species within an extremely small distribution range, on the right bank of the Volga River, most likely representing a contemporary refugium. The genealogy of plastid markers in P. volgarica suggests that this species is ancestral to P. eversmanniana, another narrow endemic from the S Urals.


Introduction
The genus Potentilla L. comprises more than 300 species distributed worldwide in temperate areas and in mountainous regions in the tropics [1][2][3]. Hybridization, polyploidy, and apomixis are not rare among its species [4][5][6][7][8][9][10], which make the taxonomy of the group very complicated. The existing phylogenies of Potentilla and the Potentilleae tribe are based on relatively small subsets of taxa and are still far from comprehensive [2,[11][12][13][14][15]. Even in the cases where the taxa sets of Potentilleae were quite comprehensive [2,13], some groups of Potentilla s. str. were underrepresented, including those that are the focus of our study. The taxonomy of the genus is, similarly, far from a definite assessment. Despite the existence of a number of relatively recent regional revisions and critical taxonomic accounts [7,[16][17][18][19][20][21] the taxonomic system of the genus are still based on the monograph by 1.
To assess the genetic variability of P. volgarica and P. eversmanniana and to test whether they represent two separate species.

2.
To assess the genetic distinctions of both species from P. multifida agg. sensu Soják [19] and other related species of the section Multifidae.

3.
To assess the phylogenetic relationships of P. multifida agg. species and pinpoint the origin of disjunct isolated populations of P. volgarica and P. eversmanniana in the Russian Plain and the foothills of the Southern Urals, respectively.

Plastid Data Analyses
The length of the trnH-psbA IGS varied from 376 bp to 488 bp in the ingroup (Potentilla sect. Multifidae) and from 308 bp to 416 bp in the outgroup (Potentilla species from other sections). The length of the ndhC-trnV IGS varied from 489 bp to 543 bp in the ingroup and from 506 bp to 585 bp in the outgroup. The length of the concatenated alignment was 1215 bp. The length of the alignment after trimming with the BMGE (Block Mapping and Gathering with Entropy) v. 1.1. software [30] and manually removing three 'AT' short repeats causing homoplasy in preliminary analyses was reduced to 910 bp. The final alignment contained 56 variable positions, of which 19 positions were parsimoniously informative, 37 positions were autapomorphic, and 64 sites were alignment gaps treated as missing data in further analyses.
This trimmed alignment was used for a maximal likelihood analysis (Figure 1), which resulted in the inclusion of sequences of P. multifida, P. tergemina, P. anachoretica, P. arctica Lehm., P. agrimonioides, P. aphanes, P. jenissejensis, P. ornithopoda, P. approximata Bunge, P. verticillaris Stephan ex Willd., P. volgarica, P. eversmanniana, and P. nivea L. into the ingroup. The latter species initially was regarded as an outgroup member. The tree was poorly resolved and sequences of many samples were placed on short or zero length branches. That was indicative of incomplete lineage sorting within and among the taxa of Potentilla analyzed and justified the use of haplotype networks to reconstruct the relationships among them [31,32].
At the second stage, we reduced the alignment by excluding all the distantly related sequences of the outgroup and analyzed the genealogical relations of the ingroup sequences with the statistical parsimony approach realized in the TCS software. Indels were treated as missing data. The program calculated the 95% parsimony limit of 12 mutational steps and collapsed the sequences into 40 haplotypes united into a single network. Twelve of them were placed into the network as internal haplotypes, connected to two or more neighboring haplotypes. Accordingly, 28 haplotypes were tip haplotypes, connected to a single neighboring haplotype [33]. Further twenty-eight haplotypes were calculated by the program and included into the network as missing hypothetical haplotypes. The network has no loops and is shown in Figure 2. It comprises 12 internal haplotypes designated with capital letters A to L, and 28 tip haplotypes designated with letters and figures. Of course, these haplotype names are conventional and their choice is, to a large extent, driven by the convenience of further network interpretation. The network is unrooted and can be considered as a combination of several variously related haplotype groups, each encompassing internal haplotypes and tip haplotypes derived from them. Their geographical distribution and correspondence to morphological taxa is shown in Figure 3. Terminal names within the ingroup are followed by corresponding haplotype designations. Bootstrap support higher than 50% is indicated above branches. Different haplotype groups are highlighted with colors corresponding to those in Figure 2. Bootstrap support higher than 50% is indicated above branches. Different haplotype groups are highlighted with colors corresponding to those in Figure 2.
The second group represents haplotypes descending from internal haplotype G, differing from haplotype A by one mutational step. They include tip haplotypes G1 to G10 differing from haplotype G by one to four mutational steps, and a lineage comprising internal haplotype H and tip haplotype H1. Internal haplotype G occurs among samples of P. agrimonioides from the Caucasus, P. multifida from Southern Siberia, and P. tergemina from different parts of its area, including the samples occurring as weeds carried along railways in the European part of Russia and Ukraine. As in the previous case, tip haplotypes are rare and found in solitary samples of several species. Haplotypes G1, G5, G6, G8, G9 are found in P. tergemina, haplotypes G2, G3, and G4 occur in P. multifida, and haplotype G7 is characteristic of one sample of P. agrimonioides from the Northern Caucasus. Haplotypes of the H-H1 lineage are exclusively characteristic of P. anachoretica from Wrangel Island off the coast of Chukotka in the Arctic Ocean. The first group of haplotypes comprises tip haplotypes related to internal haplotype A. Most of them differ from internal haplotype A by a single mutational step (A2-A6). Haplotype A1 differs from A by two mutational steps and haplotype A7 differs by six mutational steps. Haplotype A is also ancestral to four lineages encompassing both internal and tip haplotypes, B-B1, C-C1, D-D1, and E-F-F1. The internal haplotype A occurs in populations of P. volgarica ( Figure 3) but was also found in single samples of P. agrimonioides from the Altai Mountains., P. jenissejensis from Tyva, and P. nivea from the Caucasus. As to tip haplotypes and the clades descendant from the internal haplotype A, the pattern is more complex. Haplotypes A1-A5 were found in single samples of P. multifida from Tyva, P. aphanes from Southern Tajikistan, P. volgarica from Saratov Province, P. anachoretica from the Taimyr Peninsula, and P. nivea from the Western Caucasus, respectively. Haplotype A6 was uniquely found in multiple samples of P. arctica from the shore and islands of the White Sea. Haplotype A7 was found in P. verticillaris from the shore of Lake Baikal. The lineages B-B1, C-C1, D-D1, and E-F-F1 appeared to be species specific for P. anachoretica, P. jenisseiensis, P. agrimonioides, and P. ornithopoda, respectively. The second group represents haplotypes descending from internal haplotype G, differing from haplotype A by one mutational step. They include tip haplotypes G1 to G10 differing from haplotype G by one to four mutational steps, and a lineage comprising internal haplotype H and tip haplotype H1. Internal haplotype G occurs among samples of P. agrimonioides from the Caucasus, P. multifida from Southern Siberia, and P. tergemina from different parts of its area, including the samples occurring as weeds carried along railways in the European part of Russia and Ukraine. As in the previous case, tip haplotypes are rare and found in solitary samples of several species. Haplotypes G1, G5, G6, G8, G9 are found in P. tergemina, haplotypes G2, G3, and G4 occur in P. multifida, and haplotype G7 is characteristic of one sample of P. agrimonioides from the Northern Caucasus. Haplotypes of the H-H1 lineage are exclusively characteristic of P. anachoretica from Wrangel Island off the coast of Chukotka in the Arctic Ocean.
The first group of haplotypes comprises tip haplotypes related to internal haplotype A. Most of them differ from internal haplotype A by a single mutational step (A2-A6). Haplotype A1 differs from A by two mutational steps and haplotype A7 differs by six mutational steps. Haplotype A is also ancestral to four lineages encompassing both internal and tip haplotypes, B-B1, C-C1, D-D1, and E-F-F1. The internal haplotype A occurs in populations of P. volgarica ( Figure 3) but was also found in single samples of P. agrimonioides from the Altai Mountains., P. jenissejensis from Tyva, and P. nivea Internal haplotypes I and J are three and four mutational steps from haplotype G, respectively, and are found in two samples of P. tergemina, from a railway in Kiev (Ukraine) and a roadside plant in the Southeastern Urals, respectively.
Internal haplotype K is distanced from haplotype J by four mutational steps. It is central for the third group of haplotypes represented by tip haplotypes K1-K4 and a clade L-L1-L2. Haplotype K is characteristic of P. eversmanniana from the Southern Urals. However, it is also found in a sample of P. multifida from the Altai Mountains, a sample of P. anachoretica from Wrangel Island, and in four samples of P. volgarica from a single locality (a chalk hill near Novaya Yablonka) in Saratov Province. Tip haplotypes were found in solitary samples of P. volgarica (K1), P. eversmanniana (K2), P. multifida (K3), and P. agrimonioides (K4). The haplotypes of the lineage L-L1-L2 were exclusively found among samples of P. volgarica.
To root the network, we reduced the alignment, deleting all identical sequences among accessions of the same species, and analyzed it together with outgroup sequences with the maximal likelihood approach in raxmlGUI to reconstruct the species tree. The resulting tree is shown in Figure 1. Separate species are represented here by one to several accessions corresponding to the haplotypes revealed with TCS. Generally, the tree is congruent with the statistical parsimony-based haplotype network, yet the basal node is weakly supported and unresolved. The basal node forms a polytomy in which accessions from P. volgarica, P. nivea, P. jenissejensis, and P. agrimonioides corresponding to A haplotype sequences are positioned on zero (or very close to zero) length branches, thus matching the internal position of the A haplotype in the network ( Figure 2). Eight more terminals of the basal polytomy are positioned on non-zero length branches, representing accessions of P. agrimonioides (from the Altai Mountains), P. anachoretica, P. aphanes, P. arctica, P. multifida, and P. volgarica. The basal polytomy also contains three highly (98%) to moderately (81%) supported clades corresponding to the three haplotype lineages derived from the internal haplotype A in the network. These are the clades of P. anachoretica (B and B1 haplotypes), P. jenissejensis (C and C1), and P. ornithopoda (E, F, and F1).
One of the lineages present in the network is not supported by the tree (D-D1). The major clade derived from the basal polytomy is weakly supported and includes all the remaining samples. It also forms a polytomy, and encompasses three samples of P. agrimonioides, P. multifida, and P. tergemina corresponding to the internal haplotype G and positioned on zero length branches. Nine more terminals emerge from the polytomy on mostly short branches representing accessions of P. multifida, P. tergemina, P. approximata, and P. agrimonioides (the Caucasus). In addition to these, the polytomy contains two moderately supported clades, one encompassing two accessions of P. anachoretica and P. agrimonioides (86%) and the other uniting the remaining samples (79%). The latter clade includes a basal grade of two accessions of P. tergemina (representing internal haplotypes I and J of the network) and another polytomy uniting samples of P. eversmanniana, P. anachoretica, P. volgarica, P. multifida, and P. agrimonioides sharing internal haplotype K and all its derivatives.

Nuclear ITS Data Analyses
We managed to sequence the nuclear ribosomal ITS region from only a subset of samples sequenced for plastid IGS regions (Appendix A, Table A1). Readable parts of the ITS region varied in length from 390 to 521 bp. The alignment length was 529 bp, starting from the motif TTGTCGAA to the motif GAGGCT(T/-)CC, without any major gaps. Thirty-six sequences of the ingroup and five of the outgroup had 1-11 positions with ambiguities due to double peaks in electrophoregrams indicating probable heterozygosity of the samples. Altogether the alignment had 144 polymorphic sites, 71 of which had more than two variants. We did not clone sequences with ambiguities, but reconstructed possible ribotypes using the PHASE algorithm [34,35] as realized in DNAsp. Though plants under study are most probably not diploids (see Introduction), we assumed them to be diploids for the purpose of further analyses. The alignment thus obtained had two sequences per individual, representing reconstructed alleles or ribotypes. We analyzed it using the ML approach in raxmlGUI. The resulting best tree was not fully resolved and many terminal branches were of zero length (Appendix A, Figure A1). The tree was converted to cladogram format for convenience of interpretation ( Figure 4). Two different alleles of the same accession are designated with Figures 1 and 2 after a hyphen character in a terminal name. The ingroup forms a monophyletic clade with 100% bootstrap support, with both accessions of P. nivea included into the ingroup. Though deeper nodes of the tree are mostly unsupported, it is notable that conspecific accessions here form monophyletic groups with few exceptions. The basal grade includes both alleles of the first accession of P. nivea from the Caucasus and the first alleles of three accessions of P. anachoretica from Wrangel Island. Their counterparts with all the remaining accessions of this species from Wrangel Island constitute clade I. Clade II comprises two accessions of P. agrimonioides from the Caucasus representing both alleles of accession AGR2, and one of the alleles of accession AGR9. The second allele of this accession appears to be in the next (not numbered) clade of the grade together with an allele of the second accession of P. nivea. Clade III unites all the accessions of P. volgarica. This clade is weakly supported (53%), yet most of its terminal subclades have moderate to high support. Clade IV unites most of the accessions of P. agrimonioides, both from the Caucasus and the Altai Mountains, and all the three accessions of P. jenissejensis, which form a separate subclade. Clade V unites accessions of several species and includes three major subclades. The subclade Va (unsupported) unites alleles of three accessions of P. eversmanniana; the subclade Vb (69% support) unites three accessions of P. multifida and one accession of P. anachoretica (both alleles); the subclade Vc unites most of the accessions of P. arctica from the Kola Peninsula and one of the alleles of P. aphanes from Tadjikistan. Clade VI is the most heterogenous and includes accessions of P. tergemina, P. arctica, P. multifida, P. ornithopoda, and two alleles of P. aphanes and P. nivea with their counterparts in different clades.  Table A1 and Figures 1 or 2 indicating reconstructed alleles. Major clades corresponding to species are designated with Roman numerals. Bootstrap support higher than 50% is indicated above branches.  Table A1 and Figure 1 or Figure 2 indicating reconstructed alleles. Major clades corresponding to species are designated with Roman numerals. Bootstrap support higher than 50% is indicated above branches.

Potentilla volgarica and P. eversmanniana Population Structure Analyses
We specifically analyzed populations of P. volgarica and P. eversmanniana to assess intra-and interpopulation genetic variability. Populations of P. volgarica appear to be extremely polymorphic in terms of plastid data and slightly less so in terms of nuclear ITS data, while populations of P. eversmanniana show low to zero variability in haplotype compositions (Table A1). The results of AMOVA analyses (Table 1) based on plastid sequences show most of the variability is between in-group populations (i.e., conspecific populations in our case). In the case of populations of P. volgarica taken separately, most of the variability is among local populations. The Mantel test demonstrated medium, but significant, correlation between genetic and geographical distances when both species are taken into consideration (r = 0.497, p = 0.000). However, correlation between genetic and geographical distances is not significant (r = 0.470, p = 0.084) when only populations of P. volgarica are considered. Populations of P. volgarica are geographically structured at a local scale, haplotypes of the most derived plastid L clade occurring at the extreme south and north of the species area, whereas haplotypes of A and K clades occupy its central part ( Figure 5).
Plants 2020, 9, x FOR PEER REVIEW

Potentilla volgarica and P. eversmanniana Population Structure Analyses
We specifically analyzed populations of P. volgarica and P. eversmanniana to assess intra-and interpopulation genetic variability. Populations of P. volgarica appear to be extremely polymorphic in terms of plastid data and slightly less so in terms of nuclear ITS data, while populations of P. eversmanniana show low to zero variability in haplotype compositions (Table A1). The results of AMOVA analyses (Table 1) based on plastid sequences show most of the variability is between in-group populations (i.e., conspecific populations in our case). In the case of populations of P. volgarica taken separately, most of the variability is among local populations. The Mantel test demonstrated medium, but significant, correlation between genetic and geographical distances when both species are taken into consideration (r = 0.497, p = 0.000). However, correlation between genetic and geographical distances is not significant (r = 0.470, p = 0.084) when only populations of P. volgarica are considered. Populations of P. volgarica are geographically structured at a local scale, haplotypes of the most derived plastid L clade occurring at the extreme south and north of the species area, whereas haplotypes of A and K clades occupy its central part ( Figure 5).

Discussion
Our results suggest P. multifida agg. comprises a number of relatively young and incompletely genetically differentiated species widely distributed in Northern Eurasia. Plastid data suggest an incomplete lineage sorting (ILS) characteristic of the group as a whole, including P. nivea, traditionally referred to as a different section Niveae (Rydb.) A.Nelson. As it is clear from the plastid species tree (Figure 1), P. nivea shares the most basal haplotype A with a number of accessions of different species of P. multifida agg. The internal basal haplotype A was abundantly sampled from populations of P. volgarica only. In addition to these, we managed to reveal the haplotype A only in two accessions of P. agrimonioides and P. jenissejensis from the Altai Mountains in Altai and Tyva Republics, respectively, and in the above-mentioned accession of P. nivea from the Northern Caucasus. At the same time, its derivative tip haplotypes and clades are widely distributed over the whole range of the group. They are however absent from the Caucasus, the Urals, and forested areas of Southern Siberia, mostly occupied by populations bearing haplotypes of derivative G and K haplotype lineages (Figures 1 and 3). The haplotype A clade members appear to prevail in more harsh environments in the north and high mountains in the south. However, generally, no clear geographic pattern can be seen in the distribution of plastid haplotype groups with several instances of distantly-related tip haplotypes occurring in the same population. The absence of a clear distribution pattern supports incomplete lineage sorting and recurrent hybridization. The picture emerging from plastid data is that, much like European P. crantzii (Crantz) Beck ex Fritsch [26], the common ancestor of P. multifida agg. occupied a wide Eurasian periglacial range during cold periods of the Pleistocene period, and contracted to the modern disjunct distribution of P. multifida agg. species with climate warming. Similarly, considerable range expansions during the Last Glacial Maximum, and corresponding range contractions during the last interglacial and Holocene periods were discovered in Sibbaldia procumbens s.l. [36] and Rosa sericea s.l. complex [37].
Quite surprisingly, nuclear ITS demonstrates a different pattern, with most conspecific accessions, notably P. anachoretica, P. volgarica, and P. agrimonioides together with P. jenissejensis, P. eversmanniana, P. multifida, P. arctica, and P. tergemina, nesting within monophyletic clades ( Figure 4). Several exceptions, where separate alleles of these species fall outside their main clades, appearing mostly in the clade VI representing all the accessions of P. tergemina, may indicate hybridization of these species with P. tergemina or its direct ancestor in the past. Notably, P. ornithopoda inferred alleles are dispersed among subclades of clade VI, while in the plastid tree its accessions form a monophyletic clade of haplotypes E, F, and F1. This may indicate a hybrid origin from unknown parents, probably unsampled in our study.
A special case represents P. volgarica. First of all, it is extremely genetically diverse as to plastid haplotypes. Most of the other species considered in this study possess few usually closely related haplotypes, or even are monomorphic, as P. arctica (A6), or nearly monomorphic, as P. eversmanniana (K, K2). P. volgarica is the only species represented by seven different haplotypes, A, A3, K, K1, L, L1, and L2, from two distant haplotype groups (A and K). Though just a subset of samples was sequenced for nrITS, all but one accession appeared to be homozygous, a situation reversed in other species where heterozygotes predominate. This also refers to its probably closest relative, P. eversmanniana from the Southwestern Urals. We sampled two local populations of this species (Table 1) and revealed them nearly monomorphic as to plastid haplotypes-all the plants possessed haplotype K, while its close derivative tip haplotype K2 was found in a singe plant. Contrary to that, all three specimens sequenced for ITS appeared to be heterozygotous. Moreover, populations of P. volgarica are mostly represented by homozygotes. We have studied seven local populations of P. volgarica from all localities of this species known so far ( Figure 5). This mosaic pattern of plastid haplotype diversity together with predominant homozygosity of populations by ITS and results of AMOVA and Mantel test is suggestive of facultative apomixis in this species.
The unusual pattern of genetic diversity in P. volgarica may be explained by long persistence of the species in its current, extremely small distribution area, which probably behaves as a contemporary refugium. Ecologically, this species is restricted to steppe on hills with chalk outcrops. Potentilla volgarica populations, especially those in the central part of the area, deserve protection because they harbor most of the species' genetic polymorphism. Potentilla eversmanniana needs further study with larger sampling. However, the low genetic diversity observed in the present study suggests that it may be more vulnerable to habitat disturbance and climate change, than P. volgarica. Moreover, the plastid data suggest that it diverged from P. volgarica in the Pleistocene period, when this species probably had a wider distribution area.

Taxon Sampling
The plant material used in the present study was sampled from five local populations in the Saratovskaya Province of Russia (P. volgarica) and from two local populations in the Republic of Bashkortostan of Russia (P. ewersmanniana) in May-June 2019. Additional samples of both species, as well as P. arctica (a two population series in MHA and MW), P. anachoretica (a population series in MW and MHA), P. agrimonioides, P. aphanes Soják, P. jenissejensis (three specimens kept in MW under the name P. multicaulis), P. multifida, P. ornithopoda, P. approximata Bunge, P. tergemina (including several samples outside the natural range of the species collected at railroads in European Russia), and P. verticillaris were obtained from herbarium specimens kept at MHA, MW, and SARBG Herbaria. We determined the sampled plants using the keys to species of Potentilla in "Flora Europae Orientalis" [7], "Monographie der Gattung Potentilla" [1], J. Soják's critical papers [18][19][20], and V. Kurbatsky's [38] keys to the species of Potentilla of Asian Russia based on morphological characters. Two samples of P. eversmanniana from MW collected in Sverdlovsk Province (the Central Urals) were redetermined as P. tergemina, and three samples of P. multicaulis Bunge as P. jenissejensis. All the samples used for the analyses are listed in Appendix A, Table A1.

Molecular Methods
Total DNA was extracted from silica-dried leaf tissue and, in some cases, from herbarium samples using the NucleoSpin Plant DNA kit (Macherey Nagel, Germany) according to the manufacturer's protocol. For the molecular phylogenetic study we used three markers: nuclear ribosomal ITS1 and two plastid intergenic spacers (IGS), ndhC-trnV and psbA-trnH. For the amplification and subsequent sequencing of ITS region, NNC-18S10 (AGGAGAAGTCGTAACAA) and C26A (GTTTCTTTTCCTCCGCT) primers were used [39]. The plastid psbA-trnH IGS was amplified using trnH (CGCGCATGGTGGATTCACAATCC) and psbA (GTTATGCATGAACGTAATGCTC) primers, and the ndhC-trnV IGS was amplified with ndhC (ATTAGAAATGYCCARAAAATATCAT) and trnV(UAC)x2 (GTCTACGGTTCGARTCCGTA) primers [40,41]. We chose these two regions after a pilot screening of several potentially variable plastid markers [40,41] with a small subset of samples. All PCR products were directly sequenced in both directions with the same primers. Primers used for PCR were synthesized and purified in PAAG by Syntol Ltd. (Moscow, Russia). Double-stranded DNA templates were obtained by polymerase chain reaction (PCR). PCR reaction mixture (20 µL) contained 1.5-2 ng of DNA, 5 pmol of each primer, 4 µL of Ready-to-Use PCR Master mix 5× MasDDTaqMIX-2025, containing a "hot-start" SmarTaq DNA polymerase (Dialat Ltd., Moscow, Russia) and 13 µL of deionized water. PCR reaction was performed on a MJ Research PTC220 DNA Engine Dyad Thermal Cycler  Table A1 (Appendix A).

Sequence Editing and Alignment
For the purposes of analysis, both plastid regions were concatenated. We were not able to sequence the ITS region for all samples (see Table A1); for many samples this region was sequenced incompletely. Sequences were aligned using MAFFT [42,43] with the accurate strategy L-INS-i and modified manually using BioEdit 7.0 [44]. ITS and concatenated chloroplast regions were analyzed separately. Since the plastid alignment had multiple indels, some of which being dubiously aligned, we used the BMGE (Block Mapping and gathering with Entropy) v. 1.1. software [30] with the default -t option to trim the alignment to remove gaps and phylogenetically uninformative and homoplasious positions.

Data Analyses
Phylogenetic reconstruction was performed in RAxML ver. 8.2.10 using raxmlGUI 2.0 beta [45][46][47][48]. Bootstrap values are based on 100 replicates (Fast bootstrap option). The program searched for trees with the maximum likelihood approach under the GTRGAMMA model with parameters calculated by the program. Phylogenetic relationships among the cpDNA haplotypes were reconstructed using statistical parsimony analysis as implemented in TCS v1.2 [49] with alignment gaps treated as missing data. We edited the resulting trees in FigTree v.1.4.3 [50] and finally elaborated all the figures in InkScape v.0.48.2 (https://inkscape.org).
To root the trees we downloaded from GenBank eight complete chloroplast genomes of different species of Potentilla and used only two regions from these cp genomes, ndhC-trnV and psbA-trnH. Additionally, we sequenced two specimens of P. nivea as a possible close outgroup [14].
The programs DNAsp v.6 [51] Figure A1. Maximum likelihood tree of Potentilla multifida agg. based on ITS data. Terminal names within the ingroup are followed by accession designations as in Table A1 and Figures 1 or 2 indicating reconstructed alleles. Lower-case letters after P. volgarica terminals indicate ribotypes of this species. Different species are highlighted with different colors. Figure A1. Maximum likelihood tree of Potentilla multifida agg. based on ITS data. Terminal names within the ingroup are followed by accession designations as in Table A1 and Figure 1 or Figure 2 indicating reconstructed alleles. Lower-case letters after P. volgarica terminals indicate ribotypes of this species. Different species are highlighted with different colors.