Phylogeny and biogeography of the Pleistocene Holarctic steppe and semi-desert goosefoot plant Krascheninnikovia ceratoides

Krascheninnikovia ceratoides (Chenopodiaceae) is a steppe and semi-desert plant with two subspecies, K. cer- atoides subsp. ceratoides , which is widespread in Eurasia, and K. ceratoides subsp. lanata , which grows in western and central North America. A few disjunct populations of K. ceratoides subsp. ceratoides are found in Anatolia, Europe and North Africa to the west of its otherwise continuous Eurasian distribution. To understand the evolutionary history of this characteristic steppe and semi-desert plant, we analysed its phylogeny and biogeo- graphy. We sequenced several loci including ITS, ETS and the chloroplast intergenic spacer regions atpB-rbcL , rpl32-trnL and trnL-trnF to establish a time-calibrated phylogeny and reconstruct intraspeci ﬁ c relationships. Furthermore, we identi ﬁ ed the ploidy level of individuals. While diploid, tetraploid and hexaploid individuals have been reported in the literature, we were only able to ﬁ nd diploids and tetraploids. The diploids were found in the east of Mongolia, Kazakhstan, Kyrgyzstan, Russia and the USA. The tetraploids were located in the west of Mongolia, Kazakhstan, Tajikistan, Russia, and Europe. Populations were uniformly di- or tetraploid. Our results indicate that the species may have spread from the area of Mongolia, northern China and Middle Asia (i.e., the Altai Mountains region) in two opposite directions – on the one hand, diploids migrated to the east, to eastern Asia and North America, and on the other hand diploids and tetraploids migrated to the west, to western Asia and Europe. Fossil-calibrated gene trees were used to estimate the age of the species. Diversi ﬁ cation within the species is probably of Pleistocene age. Our dated analysis indicates that the ﬁ rst split among extant lineages of the species took place in the Early Pleistocene (Gelasian). The spread of the main lineages is likely related to major phases of steppe and semi-desert expansions during glacial periods of the Pleistocene.


Introduction
The goosefoot genus Krascheninnikovia Gueldenst. (Chenopodiaceae) is represented by a single species, Krascheninnikovia ceratoides (L.) Gueldenst., that inhabits shortgrass steppes, desert steppes and semi-deserts of the Holarctic temperate region with extensions into the subtropical region throughout the Northern Hemisphere (Heklau and Röser, 2008;Heklau and von Wehrden, 2011;Pfadenhauer and Klötzli, 2014). Its main distribution broadly corresponds to the Irano-Turanian floristic region (Takhtajan, 1986), but its entire range extends far beyond. The wide distribution of the species and high levels of morphological variation have resulted in a confusing and contradictory taxonomy (Rechinger, 1963;Komarov, 1964;Täckholm, 1974;Welsh et al., 1987;Davis, 1988;Tutin et al., 1993;Wu et al., 1994;Heklau, 2006;Holmgren, 2008). However, Heklau and Röser (2008) recognised two subspecies within it based on both morphological and genetic data: the Old World (mainly Eurasian) subspecies ceratoides and the New World subspecies lanata (Pursh) Heklau. Krascheninnikovia ceratoides subsp. ceratoides has a broad and continuous distribution in Asia from Mongolia and Inner Mongolia (China) in the east to the southern Volga and Don region (Russia) in the west (Fig. 1). To the south, it reaches the Himalayas and Iran. To the west of this continuous distribution, small and isolated exclaves exist in Anatolia (Turkey), Europe (Austria, Romania, Spain and Ukraine) and North Africa (Egypt and Morocco). Krascheninnikovia ceratoides subsp. lanata is distributed in central and western Canada, USA and Mexico ( Fig. 1).
It is assumed that the formation of the Eurasian semi-deserts and steppes started in the Neogene in Asia (Willis and McElwain, 2002). During the Quaternary, the semi-deserts and steppes faced extensive range shifts, contractions and expansions (Frenzel, 1968;Lang, 1994;Franzke et al., 2004;Friesen et al., 2016). During periods of cooling and the peaks of the Pleistocene glaciations, semi-desert and steppe communities reached their maximum extent. The sea level dropped, and the climate was very arid (MacNeil et al., 1943;Frenzel, 1968;Pfadenhauer and Klötzli, 2014). In contrast, during the warm and humid Pleistocene interglacial periods, the forests recolonised the dry and open areas (González-Sampériz et al., 2010). Thus, the distribution of steppe and semi-desert plant species could correlate with the expansions and contractions of the steppe and semi-desert biomes during the Pleistocene.
The disjunct distribution pattern exhibited by Krascheninnikovia ceratoides subsp. ceratoides has intrigued researchers for a long time (Willkomm and Lange, 1870;Braun-Blanquet and Bolós, 1957;Costa Tenorio et al., 2000). Three main explanations have been proposed. The first hypothesis suggests that the extant isolated populations of K. ceratoides subsp. ceratoides in Europe and North Africa are relicts of a pre-Quaternary (Miocene or Pliocene) stepping-stone migration of the plant from Asia to the west (Bolós, 1951;Braun-Blanquet and Bolós, 1957;Sainz Ollero et al., 1996;Costa Tenorio et al., 2000), possibly facilitated by the desiccation of the Mediterranean basin during the Messinian salinity crisis (6.0 Mya; Manzi et al., 2013). The second hypothesis proposes that the disjunct populations are the result of recent migration and/or long-distance dispersal in the Pleistocene, in keeping with the explanation proposed for vicariant distributions in other plants and animals (e.g., insects, Ribera and Blasco-Zumeta, 1998; Microcnemum coralloides, Kadereit and Yaprak, 2008; Smilax aspera, Chen et al., 2014;Cheirolophus, Vitales et al., 2014). Further support for this hypothesis comes from recent phylogeographic studies, which suggest that some southern European plants that were previously considered 'Tertiary' relicts might actually be the result of recent (Quaternary) dispersal/ speciation events (e.g., Microcnemum coralloides, Kadereit and Yaprak, 2008;Krascheninnikovia ceratoides in Spain, Pérez-Collazos and Catalán, 2007;Ferula loscosii, Pérez-Collazos et al., 2009;Ruppiaceae, Triest and Sierens, 2014). The third hypothesis proposes a more recent introduction of the plant into the western Mediterranean as fuel and animal forage, probably during the Arab invasions in the eighth century (Willkomm, 1896;Costa Tenorio et al., 2000). Notably, all three hypotheses point to Asia as the centre of origin of the taxon based on the large number of populations and wider distribution of K. ceratoides subsp. ceratoides in that area. The colonisation of North America from Asia by K. ceratoides across the Bering Strait has also been suggested (Heklau and Röser, 2008), but no biogeographic analysis has tested this hypothesis yet.
In this study we used genome size measurements as well as DNA sequencing of nuclear and chloroplast regions to reconstruct the phylogeny and the biogeography of Krascheninnikovia. Specifically, we aimed to (1) determine the ancestral area of the genus, (2) establish divergence times to test hypotheses to explain its disjunct distribution; specifically whether it is the result of a Miocene/Pliocene (Bolós, 1951;Braun-Blanquet and Bolós, 1957;Costa Tenorio et al., 2000) or Pleistocene (Pérez-Collazos and Catalán, 2007) migration or Holocene introduction by man (Costa Tenorio et al., 2000), and (3) evaluate the origins of the polyploid populations and the potential recurrence of polyploidisation in the two Holarctic continents.

Plant material
Leaf material of Krascheninnikovia ceratoides was collected in the field or taken from herbarium material to represent the entire range of distribution of the species (Fig. 1, Appendix A). In total, 180 accessions from 153 different populations of K. ceratoides were studied with the nuclear ITS and 73 accessions with the nuclear external transcribed spacer (ETS). A reduced pool of samples (92 accessions representing 92 populations) was used for the chloroplast analysis. For phylogenetic analysis of the ITS, ETS and chloroplast data, we included samples of Axyris and Ceratocarpus from the National Center for Biotechnology Information (NCBI) GenBank and field collections as outgroups. For dating analyses, sequences of 113 additional taxa from GenBank were included in order to place Krascheninnikovia in a broader phylogenetic context. Those taxa were chosen based on Chenopodiaceae phylogenies of Heklau and Röser (2008), Kadereit et al. (2003Kadereit et al. ( , 2010 and Di Vincenzo et al. (2018).

Genome size measurements
The genome size and ploidy level of 108 samples of Krascheninnikovia ceratoides were measured by flow cytometry using the CyStain PI Absolute P kit (Sysmex, Görlitz, Germany) according to the manufacturer's manual with leaf tissue of the sample and an internal standard: Pisum sativum L. 'Kleine Rheinländerin' (2C = 8.8 pg; Doležel and Greilhuber, 2010) or Petroselinum crispum (Mill.) Fuss (2C = 4.46 pg; Yokoya et al., 2000). The nuclear DNA content was measured in a CyFlow Space (532 nm diode laser; Sysmex) until at least 5000 measurements were collected.
PCR amplification conditions for each genomic region are indicated in Appendix B. Amplification products were purified using GFX purification kits (GE Healthcare, Little Chalfont, UK) or by enzymatic treatment with Exonuclease I and FastAP Thermosensitive Alkaline Phosphatase (Thermo Fisher Scientific, Waltham, Massachusetts, USA). DNA concentration was measured using a NanoDrop spectrophotometer (NanoDrop Technologies, Wilmington, Delaware, USA) or by fluorescent quantification using the DeNovix dsDNA Broad Range Assay (DeNovix, Wilmington, Delaware, USA) according to the manufacturer's manual. Both DNA strands were sequenced. Cycle sequencing reactions consisted of 1 μL BigDye Terminator v3.1 Ready Reaction Mix (Thermo Fisher Scientific), 1.5 μL 5× sequencing buffer, 1 μL of 1 μM or 3.5 μM primer and 6.5 μL of diluted PCR product. Cycling conditions after an initial temperature of 96°C for 5/1 min were 25/35 cycles of 96°C for 10 s, 50°C for 5 s and 60°C for 4 min. Excess dye-labelled nucleotides from the sequence reaction were removed either by standard ethanol/sodium acetate precipitation followed by a run on an ABI 3730 DNA capillary sequencer (Applied Biosystems, Foster City, California, USA) at the Wolfson Wellcome sequence facility of the Natural History Museum, London, or they were cleaned up using Sephadex columns followed by a run on a 3500 Genetic Analyzer (Applied Biosystems). Sequences are available online in the European Nucleotide Archive (accession of the project: PRJEB31888; accessions of the sequences: ITS: LR537340-LR537419; ETS: LR537058-LR537120; atpB-rbcL: LR537420-LR537442; rpl32-trnL: LR585005-LR585031; trnL-trnF: LR585038-LR585064).

Data analysis
Forward and reverse sequences were assembled and edited using Sequencher ver. 4.2.2 (Gene Codes Corporation, Ann Arbor, Michigan, USA) or Geneious ver. 11.1.5 (Biomatters, Auckland, New Zealand). Consensus sequences were aligned using PASTA for ETS and the atpB-rbcL intergenic spacer (Mirarab et al., 2015) or MAFFT using the Q-INSi strategy for ITS (Katoh et al., 2017). ITS, ETS and atpB-rbcL on the one hand and rpl32-trnL and trnL-trnF on the other hand were sequenced for different, mostly non-overlapping sets of samples. To test for data combinability, therefore, a partition homogeneity test was conducted in PAUP* ver. 4.0beta10 (Swofford, 2002) on the following sets of data: (1) rpl32-trnL and trnL-trnF; (2) ITS1, ITS2 and ETS; and (3) atpB-rbcL, ITS and ETS. Partitioned and combined maximum likelihood (ML) analyses were performed using the software IQTREE . The model of nucleotide substitution was selected for each dataset using jModelTest ver. 2.1.10 (Darriba et al., 2012) based on the BIC value. The best fit model was SYM + G for ETS, ITS1 and ITS2, SYM + G+I for atpB-rbcL (but SYM + G was selected for simplicity) and GTR + G for the other chloroplast regions. IQTREE computed 100 random Maximum Parsimony (MP) trees, searching for best-scoring ML trees and estimating branch support for the best tree from 1,000 bootstrap replicates using the ultrafast bootstrap option. The distance matrix produced by the IQTREE analysis was visualized using SplitsTree (Huson and Bryant, 2006).
TCS (Clement et al., 2000) was used to establish a haplotype network for the combined chloroplast rpl32-trnL and trnL-trnF matrix within a statistical parsimony framework. The analysis was performed with a 95% parsimony connection limit.

Divergence time estimation
We used only the ITS1, ITS2 and ETS regions for divergence time estimation given the low variability of the 5.8S rDNA and the atpB-rbcL intergenic spacer which would lead to high variance in age estimates at lower taxonomic levels. Molecular dating analyses of the ITS1-ITS2-ETS matrix were performed using BEAST ver. 2.5.2 (Bouckaert et al., 2014). Each of these regions gave similar age estimates when used alone. In order to implement a primary calibration to date the times of divergence among Krascheninnikovia ceratoides samples, we included 113 taxa of the Chenopodioideae clade, where the genus is nested in (Kadereit et al., 2005). The Paleocene pollen fossil Chenopodipollis multiplex dated at 65-56 Mya was used for calibrating the crown group of the Chenopodioideae (Kadereit et al., , 2005(Kadereit et al., , 2010. We used a uniform prior to constrain the age of the root node to the age of this fossil. The interface BEAUti was used to create the BEAST input file using the following parameters: (1) symmetric models with unequal rates but equal base frequencies (SYM; Zharkikh, 1994) with gamma distributed rate variation among sites with four categories of rates defined for each of the three partitions separately, (2) a relaxed molecular clock model with uncorrelated rates drawn from a log-normal distribution that was unlinked across partitions, and (3) a birth-death model of lineage diversification with random start to infer the tree topologies. The MCMC chain was run for 50 million generations, saving data every 1000 generations, producing 50,000 estimates of dates and trees. This analysis was repeated five times. Convergence statistics were analysed in Tracer ver. 1.7.1. The trees were combined using Log-Combiner ver. 2.5.2 with a burn-in of 10%. We used TreeAnnotator to produce the maximum clade credibility tree from the post-burn-in trees using median node heights and to determine the 95% probability density of ages for all nodes in the trees. The tree was visualised in FigTree ver. 1.4.3.

Ancestral range reconstruction
To define areas for the ancestral range reconstruction, we used the framework of Brummitt (2001). The framework of Takhtajan (1986) was not used, because its floristic regions are rather large in the area of interest and almost all populations would lie within two regions namely the Circumboreal and Irano-Turanian regions. Following Brummitt (2001) When identical sequences of different areas were combined, this entry was treated as present in more than one area. We used the BioGeo-BEARS package ver. 1.1.1 (Matzke, 2013(Matzke, , 2014(Matzke, , 2018 for R ver. 3.4.3 (R Core Team, 2008) to calculate the likelihoods of three commonly used models for biogeographical inference: DEC, a dispersal-extinction cladogenesis model (Ree and Smith, 2008), DIVALIKE, a likelihood version of the dispersal-vicariance analysis (Ronquist, 1997) and BAYAREALIKE, the Bayesian inference of historical biogeography for discrete areas (Landis et al., 2013). All three models were calculated with and without the parameter 'J' that permits founder-event speciation in the reconstruction. All six models (with and without J) were compared using AICc to determine the best fitting model for the dataset. We used the maximum clade credibility tree with median node heights obtained from BEAST, which we pruned so that it only included samples of Krascheninnikovia. The maximum number of ancestral areas was set to two. Given the recent diversification of Krascheninnikovia inferred from our divergence time analysis (∼2.2 Mya), no temporal scenarios were taken into account for the analysis. Two alternative dispersal hypotheses were tested: In the first hypothesis, dispersal rates were assumed equal among areas (M0), whilst in the second hypothesis, we inferred differential dispersal rates reflecting the geographic connectivity among the operational areas (M1). In this model, dispersal rates were set to 0.8 when two areas were contiguous and to 0.2 when two areas were separated by another area and between Northern America and other areas.

Ploidy
Up to five individuals per population were measured and individuals from the same population always had the same ploidy level. Our search across the Holarctic populations of Krascheninnikovia ceratoides only found diploid and tetraploid populations. Diploids have a DNA content of about 2.9 ± 0.2 pg per cell (mean ± standard deviation; N = 38), while tetraploids have about 5.6 ± 0.2 pg of DNA per cell (N = 62). The ploidy levels were corroborated by comparison of the genome size values of the population of Oberschoderlee (Austria; P028; Appendix A) retrieved in this study with chromosome counts obtained from the same, uniformly tetraploid population (Dobeš et al., 1997). Some variation in the DNA content is likely attributable to the use of dried leaf material with different storage times. In the European area, only tetraploid populations were found, while in Asia diploid and tetraploid populations coexist (Fig. 1). In the USA only four populations have been measured, which were all diploid.

Aligned DNA data matrices
The nuclear and plastid aligned data matrices comprised accessions of Krascheninnikovia and of the close genera Axyris and Ceratocarpus that were used as outgroups. The number of accessions used is indicated in Appendix C. Identical haplotypes were represented by a single accession in the phylogenetic analyses. For the dating analysis, the ITS1, ITS2 and ETS data matrices of 191 accessions with unique sequences were used -57 accessions of Krascheninnikovia, four of Axyris, five of Ceratocarpus and 125 of other Chenopodioideae genera. Of the 686 characters, 539 were variable.

Phylogenetic analyses
3.3.1. Nuclear ITS and ETS regions combined with the chloroplast atpB-rbcL region The partition homogeneity tests were not significant (p > 0.5) and therefore the tested data matrices (rpl32-trnL and trnL-trnF; ITS1, ITS2 and ETS; atpB-rbcL, ITS and ETS) were combined. Fig. 2 shows the split network of the combined matrix of ITS, ETS and atpB-rbcL obtained from IQTREE (the corresponding ML tree is shown in Appendix E). Four clades with strong support were identified. The largest clade [bootstrap support (BS) 92%] consists of predominantly tetraploids from Eurasia (populations west of the Altai Mountains). Diploid individuals from a largely overlapping area (but extending further to the east and less far west) were grouped in another clade (BS 91%), together with one tetraploid population (P069; 'Predominantly diploids of Eurasia' in Figs. 2 and 5). In the third clade (BS 75%), individuals from western Mongolia, Nepal and the Russian Altai Mountains region are grouped together. The fourth clade (BS 80%) contains individuals from eastern Mongolia, Canada and USA. This is the only clade that combines accessions of both subspecies (subsp. ceratoides and subsp. lanata). Some of the American individuals are grouped together in a highly supported clade (BS 92%), while others are mixed with Mongolian samples. Therefore, the monophyly of the subspecies is not supported. Some Asian samples [P211 (Kyrgyzstan); P209, KcMon2, KcMon4 and AM849251 (Mongolia); P242 and KcRus12 (Russian Altai Mountains); and the samples from Nepal] are not resolved in any of the well-supported clades described. The split network showed KcPak1 as the earliest diverging lineage, splitting after the outgroups Axyris and Ceratocarpus (Fig. 2). In the ML tree, KcPak1 was sister to all other accessions of Krascheninnikovia (Appendix E). Relationships within and between groups were often poorly resolved and supported. A. Seidl, et al. Flora 262 (2020) 151504 3.3.2. Statistical parsimony network of combined rpl32-trnL and trnL-trnF chloroplast regions Chloroplast haplotype network analysis showed that the same haplotype was shared by 18 individuals from Western Europe, Eastern Europe, Western Asia including Pakistan and Mongolia (Fig. 3). In total, five samples of the adjacent areas of Mongolia (KcRus6, KcRus7, KcRus8, KcRus13 and KcRus14) were included in the Eurasian and Mediterranean group. Nevertheless, seven steps differentiated most Mongolian populations from the latter group. The New World populations were resolved as a monophyletic group within which no geographical structure was evident.

Biogeography
Divergence times and ancestral range reconstructions within Krascheninnikovia are presented in Fig. 5. Based on the AICc values, we selected the DEC + J analysis using model M1 with different dispersal rates among areas. The highest percentage of likelihood indicates Mongolia and northern China and adjacent populations from southern Siberia combined with Middle Asia (DE) as the ancestral range (0.18). According to the analysis, an initial vicariance event within K. ceratoides took place shortly after ∼2.2 Mya, when the ancestor of the predominantly tetraploid Eurasian clade (in Middle Asia) separated from the ancestor of the other clades (in Mongolia and northern China). The predominantly tetraploid Eurasian lineage dispersed shortly after ∼0.7 Mya from Middle Asia to Eastern Europe and later-on also to Western Europe, Turkey and Morocco as well as to Mongolia. Approximately ∼2.0 Mya, the ancestor of the predominantly diploid Eurasian lineage separated from the other lineages that originated east of the Altai Mountains. This clade diverged within the area of Mongolia and northern China from ∼0.7 Mya onwards and its descendant lineages dispersed from this area to Middle Asia, Pakistan and Eastern Europe later-on. Two lineages are resolved in the group comprising samples from areas east of the Altai Mountains that both migrated to Northern America more recently. The first lineage dispersed to Northern America between ∼1.8 and ∼0.5 Mya. The second lineage expanded its range from the area of eastern Mongolia to Northern American between ∼1.4 and ∼0.7 Mya. Thereafter, populations in  KcSp1, KcSp2, KcSp4, KcSp5, KcSp6, KcAus1, KcAus3, KcRus1, KcRus2, KcRus3, KcRus4, KcRus6, KcRus7, KcRus13, KcRus15, KcRus16, KcEgy2, KcTur2, KcTur3, KcTur5. A. Seidl, et al. Flora 262 (2020) 151504 Northern America and in eastern Mongolia became isolated from each other. From the Mongolian and northern Chinese ancestral range of the most recently evolved group, new dispersals to Middle Asia, Pakistan and Nepal were also inferred.

Origin and age of Krascheninnikovia
To date, only one study had focused on the phylogenetic relationships, morphological variation and taxonomy within Krascheninnikovia, but no attempt to estimate the timing of the origin of the genus was performed (Heklau and Röser, 2008). Our focus here is on the detailed biogeography and evolutionary history of Krascheninnikovia. To this aim, we reconstructed the phylogeny and the biogeography of Krascheninnikovia and placed it in a temporal context by molecular dating.
Divergence times estimated in this study are consistent with other studies. Our ITS1-ITS2-ETS phylogenetic tree (Appendix F) showed a similar topology to that presented in Kadereit et al. (2003Kadereit et al. ( , 2005, and our BEAST analysis dated the split between Krascheninnikovia/Ceratocarpus and Axyris at 21.8 ± 6.2 Mya. A Chenopodioideae phylogeny using a relaxed molecular clock implemented in BEAST based on atpB-rbcL sequences provided a split node age of 21.2 Mya between K. ceratoides and Axyris prostrata (Kadereit et al., 2010). A BEAST analysis based on ITS sequences by Di Vincenzo et al. (2018)  Our data support the monophyly of Krascheninnikovia, Ceratocarpus and Axyris as proposed by other authors (Heklau and Röser, 2008;Kadereit et al., 2010; Appendix F) as well as the sister group relationship of Krascheninnikovia and Ceratocarpus. Kadereit et al. (2010) supported this relationship based on flower and fruit morphology given that the female flowers of Krascheninnikovia and Ceratocarpus lack a perianth and are enclosed by two bracts that are persistent in fruit, producing only one type of fruit/seed, whereas Axyris presents a threeparted simple perianth and exhibits heterocarpy and heterospermy (producing seeds differing in the thickness of the testa). This is at odds with the conclusions of Heklau and Röser (2008) who proposed a closer relationship of Krascheninnikovia and Axyris based on a detailed morphological, morphometric and molecular analysis.

Genetic relationships and diversity within Krascheninnikovia
The chloroplast and nuclear regions tested indicate a high genetic homogeneity among populations of Krascheninnikovia ceratoides. Even with the four polymorphic regions used (ITS, ETS, rpl32-trnL, trnL-trnF) low genetic diversity was detected, and several haplotypes were shared among wide areas of distribution (Fig. 3). Low levels of genetic variability in the genus were reported by Heklau and Röser (2008) in a study with a smaller number of populations. They did not find support for the differentiation of species usually accepted for Eurasia (K. ceratoides), Asia (K. arborescens, K. eversmanniana and K. lenensis) and North America (K. lanata). Based on this evidence, they proposed that K. arborescens, K. eversmanniana, K. lenensis and K. pungens should be considered conspecific with K. ceratoides (Heklau and Röser, 2008). Our phylogeny shows the nested position of K. lanata within K. ceratoides (Fig. 5), supporting its subordination to K. ceratoides sensu stricto. Whilst ITS showed low genetic diversity in this study and in the study of Heklau and Röser (2008), hypervariable molecular markers such as ISSRs have detected higher levels of genetic diversity in populations of Krascheninnikovia ceratoides within a small geographical area (Pérez-Collazos and Catalán, 2007). These higher levels of genetic diversity might explain the enormous morphological variation detected within the species (Heklau and Röser, 2008;Heklau and von Wehrden, 2011).
Although five different ploidy levels, namely diploid, triploid, tetraploid, pentaploid and hexaploid, and the existence of populations with various ploidy levels have been mentioned in the literature, we only found purely diploid or tetraploid populations. However, it should be noted that the mixed population reported in the literature was found in Kyrgyzstan (Rubtsov et al., 1989) and our samples from this area could not be measured due to their age.

Evolution and diversification of Krascheninnikovia
Central Asia has been inferred as the centre of diversification of Krascheninnikovia based on (1) the elevated number of populations and the abundance of individuals (Braun-Blanquet and Bolós, 1957), (2) the high and variable ploidy levels, where di-, tetra-, and hexaploid individuals co-occur in the same population (Rubtsov et al., 1989; in contrast to the exclusively tetraploid Western European populations or the diploid and tetraploid Northern American populations (Castroviejo and Soriano, 1990;Sainz Ollero et al., 1996;Domínguez et al., 2001;Holmgren, 2008;Li, 2008), and (3) the elevated levels of morphological variation of the central Asian populations, which in several cases have been cause of mistakes in identification and the description of new species, leading to a complicated taxonomy and a large number of synonyms of K. ceratoides (Rechinger, 1963;Komarov, 1964;Täckholm, 1974;Welsh et al., 1987;Davis, 1988;Castroviejo and Soriano, 1990;Tutin et al., 1993;Wu et al., 1994;Heklau, 2006;Holmgren, 2008). Our ancestral range analysis suggested Mongolia and northern China in combination with Middle Asia, as the ancestral area for the genus. The ancestral range with the second highest support was Mongolia and northern China in combination with northern America. Since the next relatives of Krascheninnikovia, Ceratocarpus and Axyris are distributed in Asia, it seems plausible that Krascheninnikovia originated in the area of Mongolia, northern China and Middle Asia, which agrees with our result. This region is thought to have been the cradle of Eurasian temperate grassland vegetation at the beginning of the Pleistocene (∼2 Mya; Janis, 1993;Willis and McElwain, 2002), a period corresponding to the Gelasian (Gibbard et al., 2010). Moreover, the contact zone between Mongolia, northern China and Middle Asia, the Altai Mountains, probably provided an important refugium during the Pleistocene (Fedeneva and Dergacheva, 2003;Hais et al., 2015). Interestingly, however, a sample from Pakistan (KcPak1) is resolved as the first branching lineage within the genus in the maximum likelihood tree and network ( Fig. 2 and Appendix E), but not in the Bayesian time-tree (Fig. 4, Fig. 5). This sample together with the other samples that are resolved as early branching within K. ceratoides in the ML analysis (P209, KcMon2, KcMon4 and AM849251: Mongolia; P211: Kyrgyzstan; P242 and KcRus12: Russian Altai Mountains; KcNep1 and KcNep2: Fossil-based calibration imposed to the Chenopodiaceae Paleocene pollen fossil of Chenopodipollis multiplex was used to calibrate the root node at 65-56 Mya, following Kadereit et al. (2003Kadereit et al. ( , 2005Kadereit et al. ( , 2010 and Di Vincenzo et al. (2018;see text). The chronogram shows the divergence time estimates and 95% confidence intervals for each node. Numbers are expressed in million years (My). Identical sequences of K. ceratoides were reduced to a single sequence. The chronogram was pruned to only contain samples of Krascheninnikovia. P034_etal, P094_etal and P101_etal represent several samples with identical sequences as shown in Fig. 5.
Nepal) suggests that the origin of K. ceratoides might also have been further south like the region of the Tian Shan or Pamir Mountain ranges, also suggested for other Central Asian plant groups (e.g., Lagochilus; Zhang et al., 2017;Ammopiptanthus, Shi et al., 2017).
During the glacial periods of the Pleistocene, the areas of steppe expanded and allowed the expansion of the distribution of plants such as Krascheninnikovia ceratoides that are adapted to a cold and arid climate. Consequently, the radiation and diversification hypothesis of Krascheninnikovia ceratoides suggests that during the Pleistocene, when the steppe area expanded during glacial periods, populations in the contact zone between western Mongolia, China and Middle Asia migrated to (1) the west, reaching Middle Asia, Europe, Anatolia and northern Africa and (2) the east, reaching eastern Mongolia and Northern America. The first lineage consists of predominantly tetraploid Eurasian samples (PP = 0.99) that colonised Eastern and Western Europe, Anatolia and northern Africa starting about 650,000 years ago from Middle Asia. The second lineage showed a basal separation into a group of predominantly diploid Eurasian samples (PP = 0.98) and another group that originated east of the Altai Mountains (with Northern America; PP = 0.18). The group of predominantly diploid Eurasian samples colonised Siberia and Eastern Europe starting approximately 680,000 years ago from Mongolia and northern China. Interestingly, both Eurasian groups seem to have spread in Eurasia during the same glaciation events but did not mix. This could be due to the different ploidy level or because the two lineages evolved independently in different habitats during the interglacial time. The ancestral populations may have survived during the interglacial phases at various refugia in the Altai Mountains. During the glaciation events the Altai Mountains probably contained several habitats, which were ideal for steppe plants (Hais et al., 2015), and it is possible that some of these spots (like intermountain depressions) were still suitable for K. ceratoides during the interglacial periods (Fedeneva and Dergacheva, 2003).
Two migration events from Asia to America seem to have taken place, starting from Mongolian and northern Chinese ancestral populations. The biogeographic analysis indicates a dispersal (E-F in Fig. 5) and a separate range expansion (E-EF; each supported with a posterior probability of 1.00 on the time-calibrated tree; Fig. 4) from the ancestral Mongolian range across the Pacific Ocean during the Pleistocene between ∼1.8 and ∼0.5 Mya. At this time, steppe biomes were present in North America (Webb, 1977;Willis and McElwain, 2002) and the region could be reached through the Bering Strait. A Beringian land bridge was built between Siberia and North America during the arid Pleistocene glaciation phases that likely facilitated the passage of steppe plants (Briggs, 1995;Yurtsev, 2001). The ML analysis ( Fig. 2 and Appendix E), however, does not provide any bootstrap support above 50% for two separated Northern American groups. This is also consistent with the haplotype network analysis of chloroplast sequences (Fig. 3), which only showed one group of related haplotypes for the New World, where the most common haplotype is represented by Fig. 5. Biogeographic reconstruction of Krascheninnikovia ceratoides, based on the Lagrange DEC + J model. Pie charts at nodes indicate the relative probabilities for alternative ancestral ranges (codes for operational areas correspond to those indicated in Appendix D) and are marked with the node age. Inferred dispersal (x-y), vicariance (x/y), range expansion (x-xy) and peripheral isolation (xy/y) events with the highest marginal probability are shown on the tree. The maximum clade credibility tree obtained from BEAST was pruned to only contain samples of Krascheninnikovia. The matrix M1 shows the dispersal rates between the areas. The samples assigned to either of the two independent Northern American groups in the dated analysis (KL22 versus KL01 and KL19). The answer to the question of a single (Heklau and Röser, 2008) or double colonisation of Northern America thus must await the analysis of more variable genetic markers. The existence of diploids and tetraploids in both Asia and America has been interpreted as independent polyploidisation processes (Heklau and Röser, 2008).
The colonisation of the southern parts of the range in Nepal, Pakistan and Iran likely occurred from the Mongolian and northern Chinese ancestral area in several colonisation events (Fig. 5). Hexaploid individuals (2n = 54) are mainly found in the southern parts of the range of Krascheninnikovia ceratoides, namely in the Pamir Mountains of Tajikistan (Zakharjeva and Soskov, 1981), in Kyrgyzstan and China (together with diploids and tetraploids; Rubtsov et al., 1989;Yang et al., 1996) and in Iran (together with tetraploids; Ghaffari et al., 2015). Perhaps a higher (namely hexaploid) ploidy level is favourable in the dry mountain regions south of the Eurasian steppe. Possible explanations include higher drought tolerance of polyploids (e.g., Zhou et al., 2019) or, more generally, a higher tolerance of polyploids to a wider range of environmental conditions (Fawcett et al., 2009;Van de Peer et al., 2009). A similar hypothesis could also be formulated for the tetraploid European populations. This recent origin matches other studies, which show the repeated differentiation opportunities of the Pleistocene, producing several events of colonisation of new areas (dispersal) and fragmentation of ancestral areas (vicariance), due to the expansion and contraction of the populations during glacial-interglacial periods (Senecio sect. Senecio, Coleman et al., 2003;Rana, Veith et al., 2003;Dysosma vesipellis Cheng, Qiu et al., 2009; Ligularia hodgsonii Hook., Wang et al., 2013;Praomys delectorum Thomas, Bryja et al., 2014).
The colonisation of the western Mediterranean by Krascheninnikovia ceratoides has intrigued researchers for a long time (Willkomm, 1896;Bolós, 1951;Braun-Blanquet and Bolós, 1957;Costa Tenorio et al., 2000), and two natural colonisation hypotheses have been proposed (stepping-stone migration from Central Asia to the western Mediterranean in the Miocene or Pliocene versus migration or long-distance dispersal from Central Asia in the Quaternary; Braun-Blanquet and Bolós, 1957;Costa Tenorio et al., 2000;Pérez-Collazos and Catalán, 2007;Pérez-Collazos et al., 2009). Our results suggest a relatively recent colonisation of the Mediterranean Basin in the Late Pleistocene (Hurka et al., 2019). The poor branch support does not allow us to specify the number of events and the source area(s) for the colonisation of the Mediterranean with certainty. Therefore, we cannot definitely determine the type of migration. However, due to the existence of suitable (loess rich) habitats between the Altai Mountains and western Europe during the glacial periods (Lang, 1994) and habitats still suitable today, as shown by the isolated populations between Asia and the Mediterranean (Austria and Romania), a stepping stone scenario seems to be the simplest explanation. According to this hypothesis, the species expanded its distribution area to the west during the glacial periods and survived in suitable, arid habitats during the interglacial periods. The genetic differences between the Mediterranean populations could be the result of their divergence during the interglacial periods, or of alternative migration routes separated by time or place. The stepping stone scenario was suggested for other steppe plants with disjunct Mediterranean distribution areas (Kadereit and Yaprak, 2008). Nevertheless, since some plants from the Mediterranean share an identical sequence to plants of Egypt, our study could not discard the hypothesis of the eighth century introduction of K. ceratoides into the western Mediterranean as a fuel and animal forage plant brought by the Arabs (Willkomm, 1896;Costa Tenorio et al., 2000).