Himalayan Origin and Evolution of Myricaria (Tamaricaeae) in the Neogene

Background Myricaria consists of about twelve-thirteen species and occurs in Eurasian North Temperate zone, most species in the Qinghai-Tibet Plateau (QTP) and adjacent areas. Methodology/Principal Findings Twelve species of Myricaria plus two other genera Tamarix and Reaumuria in Tamaricaceae, were sampled, and four markers, ITS, rps16, psbB-psbH, and trnL-trnF were sequenced. The relaxed Bayesian molecular clock BEAST method was used to perform phylogenetic analysis and molecular dating, and Diva, S-Diva, and maximum likelihood Lagrange were used to estimate the ancestral area. The results indicated that Myricaria could be divided into four phylogenetic clades, which correspond to four sections within the genus, of them two are newly described in this paper. The crown age of Myricaria was dated to early Miocene ca. 20 Ma, at the probable early uplifting time of the Himalayas. The Himalayas were also shown as the center of origin for Myricaria from the optimization of ancestral distribution. Migration and dispersal of Myricaria were indicated to have taken place along the Asian Mountains, including the Himalayas, Kunlun, Altun, Hendukosh, Tianshan, Altai, and Caucasus etc., westward to Europe, eastward to Central China, and northward to the Mongolian Plateau. Conclusions/Significance Myricaria spatiotemporal evolution presented here, especially the Himalayan origin at early Miocene ca. 20 Ma, and then migrated westward and eastward along the Asian mountains, offers a significant evolutionary case for QTP and Central Asian biogeography.


Introduction
The Tamaricaceae contains about eighty species [1] and four genera: Tamarix, Myricaria, Reumuria, and Hololachna [2]. This family, and Frankeniaceae, are defined as the salt-gland anatomical lineage [3]. Myricaria consists of about twelve -thirteen species [4][5][6][7] and occurs in Northern Temperate zone of Eurasia, mainly along the Asian mountains. There are eight species in Himalayas, many are endemic, thus forming a center of diversity for Myricaria (see Figure 1). Desvaux (1825) established the genus Myricaria and Niedenzu [8] presented the first classification. Zhang & Zhang [4] studied Myricaria in China and recognized ten species; they presumed the Himalayas to be the center of origin based mainly on the distribution of species. Gorschkova [7] described six species belonging to two sections in the flora of the former USSR.
Another issue relevant to Myricaria systematics is the species Myricaria elegans. Ovezinrlikov & Kinzikaeva [9] erected the genus Myrtama based on this species but it caused some controversy. Zhang et al. [10] used ITS sequence data to study the relationships within Tamaricaceae and regarded Myrtama as an intermediate genus between Myricaria and Tamarix [11]. After sampling four species from Myricaria and sequencing ITS, rbcL, and tRNAs Ser (GCU) and Gly (UCC), Gaskin et al. [1] found Myrtama and Hololachna to be distinct within Tamaricaceae, as did Zhang et al. [10]. However, based on additional sequence data, Hua et al. [12] and Wang et al. [13] confirmed that Myrtama should be included in Myricaria. Sampling ten species of Myricaria and sequencing cpDNA psbA-trnH and the rpL16 intron, Liu et al. [14][15] investigated the species-level phylogeographical patterns of Myricaria in western China as well as the origin of M. laxiflora, a unique subtropical species of conservation concern from the Three Gorges of the Yangtze River in Sichuan and Hubei provinces. The Himalayas were proposed as the center of origin of Myricaria by Liu et al. [14] with the estimated age of origin 1. 46-2.30 Ma.
Closely associated with the distribution pattern of Myricaria and related taxa, the QTP and Himalayan uplift during the Neogene are hypothesized to be a major influence on organism evolution in Asia (e.g. [16][17]). Following collision of the Indian and Eurasian continents at ca. 50 Ma, the altitude and range of the QTP near the Oligocene-Miocene boundary became sufficient to trigger a reorganization of the Asian climate, as evidenced by the beginning of loess deposition in the Chinese Loess Plateau and the Junggar Basin [18][19][20]. Some evidence confirms that the central areas of the QTP were raised to present altitudes by that time [16,[21][22] and uplift of the Himalayas may have also begun at that time [22][23]. Uplift of peripheral portions of the plateau has continued at various intervals [24][25][26][27][28]. A major uplift of QTP is often suggested to have occurred at 8 Ma, which also coupled with global cooling, even though Molnar [29] considered this uplift evidence to be inconclusive. Uplift of the QTP and global climate cooling and aridification [27] have been suggested causes for the evolution of many organisms [30][31][32][33][34][35][36]. As these studies have shown, rapid diversification of lineages in the QTP resulted in the migration of some species into other temperate regions, such as Central Asia, the Arctic, the Mediterranean (Caucasus-Alps) and southern Asia. Of these, connections between the QTP and adjacent arid, more northern areas can often be discerned, for example in recent studies on Hippophae rhamnoides (Elaeagnaceae) [37], Caragana (Fabaceae) [38], and Astragalus (Fabaceae) [39]. Linkage of the QTP and Africa and/or the Mediterranean is illustrated by Begonia [40] and Uvaria (Annonaceae) [41]. An example linking the QTP and Southeast Asia is Paini (Anura: Dicroglossidae) [42].
The origin of Myricaria has been associated with the QTP and Himalayas but justification has been weak. Zhang & Zhang [4] presumed that Myricaria originated from the Himalayas, only based on species distribution of the genus, whereas same opinion conducted by Liu et al. [14] from a phylogeography. Here we attempt to examine the origin and evolution of this genus and link it to the Himalayan uplift to explain the causes of its evolutionary patterns. In addition, the classification and distribution of Myricaria are examined using molecular phylogeny and biogeography.

Taxa sampled
Twelve species (seventeen samples) of Myricaria plus seven species from the outgroups Tamarix and Reummuria served as sources of DNA material (
The protocol for amplification consisted of an initial hotstart at 95uC for 2 min, followed by 30 cycles of denaturation at 94uC for 30 s, annealing at 52uC for 30 s, extension at 72uC for 90 s, and a final extension at 72uC for 10 min. PCR products were purified using the PEG precipitation procedure [46]. These were sequenced using an ABI Prism 3770 Genetic Analyzer (Shanghai Sanggon Biological Engineering Technology & Service, Shanghai, China). Sequences were aligned using CLUSTAL X software [47], and then adjusted by hand. All gaps were treated as missing characters. Finally, a combined dataset consisting of four sequences of ITS and three cpDNA trnL-trnF, rps16, and psbB-psbH, was prepared for phylogenetic analysis.

Phylogenetic analysis and divergence time estimate
The sequence dataset from twelve species (seventeen samples) of Myricaria plus seven species of Tamarix and Reummuria yielded 3202 aligned nucleotide characters from four genes: ITS, trnL-trnF, rps16, and psbB-psbH. The incongruence length difference (ILD) test of the four gene datasets was carried out using PAUP* [48], to assess potential conflicts between the different DNA fragments. This test was implemented with 100 partition-homogeneity test replicates, using a heuristic search option with simple addition of taxa, TBR branch swapping and MaxTrees set to 1000. 0.222 of incongruence length difference (ILD) tests [48] showed that the four gene datasets were not incongruent. Estimation of phylogenetic relationships and divergence time was conducted using a Bayesian method implemented in BEAST 1.5.4, employing a relaxed clock model [49][50]. We used the uncorrelated lognormal molecular clock model with a Yule process for the speciation model, GTR+I+G for the substitution model (estimated for the dataset), and a normal distribution with SD of 1 as priors on the calibration nodes to accommodate for calibration uncertainty. Minimum ages for the two normal priors were constrained to the root of all taxa and a node respectively, the family Tamaricaceae 70 Ma, and genus Tamarix 25 Ma, with a detailed description as follows. A Markov chain Monte Carlo analysis was run for 50 million generations and sampled every 1,000 generations, and two independent runs were performed to confirm the convergence of the analysis. The stationarity of each run was examined using the effective sampling size of each parameter (.200). The last 40 million generations were used to construct the maximum clade credibility tree and the associated 95% highest posterior density distributions around the estimated node ages. Optimization of ancestral distributions Tamaricaceae root constrained. Tamaricaceae is included in the order Caryophyllales [2,51] and has no reliable macrofossil record. According to molecular dating [52][53][54], the divergence time of the order has been estimated as ca. 100 Ma. Tiffney [55] [61]. In the light of these estimates, a balanced age for Tamaricaeae could be suggested as about 70 Ma; this estimate was chosen as the family root for molecular dating.
The earliest reliable fossil record of Tamarix is Miocene, from the Yunnan province of China [62]. Most of the available fossils are from the Miocene, therefore, the genus might be hypothesized to have had an origin at least in early Miocene. However, considering its wide distribution in Europe, Africa, Asia, and North America, and our limited samples mainly from China, we conservatively assigned an age of late Oligocene-early Miocene, at ca. 25 Ma for Tamarix.

Areas
In accordance with the distribution of Myricaria species along the Asian mountains (Figure 1), we divided the distribution into six areas, namely, A: the eastern Himalayas, including the eastern QTP, the Hengduan mountains, and northern and central China; B: the western Himalayas, including the western QTP and the Pamir-Alai, Kunlun-Altun, and Hendukosh mountains; C: the Tianshan mountains and Junggar-Turan deserts; D: Altai-Siberia; E: the Mongolian Plateau; and F: Asia Minor-Caucasus-Europe. These six areas are distinct in biodiversity, vegetation, and floristics [16][17][18]30,32,[36][37][38].

Optimization of ancestral distributions
To infer biogeographical events, three methods were used: a parsimony-based procedure Diva [63], S-Diva [64] and a maximum likelihood-based DEC model (Lagrange; [65][66]). These three approaches are simultaneously considered so that to assess the relevant biogeographical processes, such as vicariance, dispersal, and extinction.
Diva. Dispersal-vicariance analysis optimizes distributions for each node of the tree by minimizing the number of assumed dispersals and extinctions, and favors vicariance events [63,67]. The Diva program reconstructs widespread ancestral distributions, restricting them to single areas. Because allopatric speciation by vicariance is the null model in Diva, vicariance and range division would always be the preferred explanation if the ancestors were widespread. To avoid inferring a widespread ancestor at the root because of the presence of widespread extant taxa, a limit of two areas was set (maxareas = 2) in Diva [63]. The phylogenetic typology of the BEAST tree ( Figure 2) was input for Diva analysis.
S-Diva. (or Bayes-Diva) [64] is a program which complements Diva and implements the methods of Nylander et al. [68] and Harris et al. [21], determining statistical support for ancestral range reconstructions using multiple trees from Bayesian analysis. This has the advantage that uncertainties in phylogenetic inference can be taken into account. One hundred Bayesian MCMC trees with the last stable typologies from BEAST, and a BEAST tree typology ( Figure 2) were input into the S-Diva program.
Lagrange. A valuable new biogeographical method is parametric likelihood analysis, with a dispersal-extinction-cladogenesis model [67], as implemented in Lagrange v. 2.0.1 [65]. This method calculates the likelihood of biogeographical routes and areas occupied by the most recent common ancestor (MRCA) for a given phylogenetic tree topology (BEAST tree, Figure 2) and the present distribution of taxa. Therefore, dispersal and vicariance of lineages, represented by connection areas, can be estimated by the probabilities. This is thus a form of MRCA area reconstruction differing from the parsimony approach of Diva.

Phylogenetic analysis and divergence time estimate
The phylogenetic tree obtained from Bayesian inference in BEAST showed that Myricaria is monophyletic and Myrtama should be included in Myricaria rather than treated as a distinct genus ( Figure 2). Within Myricaria, four clades were recognized, two corresponding to the existing sections Parallelantherae Ndz. and Renantherae Ndz, the other two represent new groups to be named as sections Alpinae and Laxiflorae (see Appendix S1). Flowers and filaments of the plants are illustrated in Figure 2, to show the characteristics of the four sections. In the present phylogenetic tree, the clades of the genus and the sections have strong support, confirming the validity of taxa at the ranks of genus and section. Section Renantherae comprises the most species in the genus, and has two subclades. The two widely distributed species of this section, M. alopecuroides and M. squamosa, are located in each subclade. The crown age of Myricaria was ca. 20 Ma, and 8.83,6.35 Ma for four sections.

Optimization of ancestral distributions
The results of the three approaches Diva, S-Diva, and Lagrange ( Figure 3) showed a consistent and strongly supported pattern, particularly at the ancestral nodes for Myricaria (AB), and among the four sections, Parallelantherae (B), Alpinae (A), Laxiflorae (A), and Renantherae with AB from Diva and S-Diva, whereas only with ABCDEF/B from Lagrange. On the whole, AB, A, and B, namely the Himalayas and the QTP, should be considered as ancestral areas in Myricaria. The events occurring in areas C, D, E, and F, were considered to be dispersals, several of which can be distinguished. M. prostrata occurs in the Himalayas, and its western Himalayan distribution was indicated to be a dispersal event from the eastern Himalayas. M. bracteata in sand lands of the Mongolian Plateau was shown to be a migrant from the eastern Himalayas, Hengduan Mountains, and Northern China. Whereas the distribution of M. germanica in Asia Minor-Caucasus-Europe was came from the western Himalayas.

Phylogenetic division of sections within Myricaria
Niedenzu [8] divided Myricaria taxa into two sections: Parallelantherae Ndz. and Renantherae Ndz. Gorschakova [7] accepted this classification system in the Flora of the USSR. Zhang & Zhang [4], however, considered that the establishment of infrageneric ranks was not appropriate due to its complicated and variable morphological characters. Therefore, in the Flora of China [5][6] there is no division of infrageneric sections.
However, our phylogenetic tree ( Figure 2) yielded a clear phylogenetic division of four sections, including two that are new. Of them the Himalayan and QTP section Alpinae, containing M. prostrata, M. rosea, and M. wardii, is characterized by the prostrate and recumbent, and with an adaptation to high altitudes of 3000-5300 m. Section Laxiflorae, comprises of only species M. laxiflora, endemic to the subtropical area of the Sichuan and Hubei provinces in eastern China. Detailed descriptions of two new sections are given in Appendix S1.
Myricaria elegans Royle was originally described from the Kunawar region of the western Himalayas [69]. Based on this species and its characters of 10 stamens, flat leaves, and no obvious style, the genus Myrtama was established [9]. Qaiser & Ali [70] named it Tamaricaria, whereas Baum [71] moved Myricaria elegans to Tamarix as T. ladachensis. As mentioned, Zhang et al. [10] and Gaskin et al. [1] accepted Myrtama at generic rank. Zhang et al. [10] considered that Myricaria elegans was an intermediate and hybrid genus between Myricaria and Tamarix, related more to Tamarix. However, Hua et al. [12] and Wang et al. [13], based on sequence data, found that it would be appropriately placed in Myricaria. Our results (Figure 1) also show that inclusion of Myricaria elegans in Myricaria is suitable, since the whole of Myricaria, including Myricaria elegans, has strong support (100%) (Figure 2). This is in accordance with evidence from morphological classification [4] and molecular phylogeny [12][13][14]. The former conclusion supporting retention of Myrtama [10] was only based on ITS sequence data, which is probably not sufficient evidence [72]. While in the phylogeny of Gaskin et al. [1] (their Figure 2), only four species were sampled, and four species as clade had a strong bootstrap support (99%) for the inclusion in Myricaria [1].

Himalayan origin, ancestral inheritance, and multidiversification in the Himalayas
Our estimated crown age of ca. 20 Ma for Myricaria (Figure 2) falls into the probable early range of the Himalayan uplift in early Miocene [22], consequently allows us to speculate that uplift of this mountain range caused the origin of Myricaria. The biogeographical analytical result of Diva, S-Diva and Lagrange, showing the combined Himalayan area AB as the ancestral area for Myricaria (see Figure 3), which also supports a Himalayan origin. The present molecular dating results are in contrast to previous phylogeographical opinion of a main divergence event at the implausible age of 1.46-2.30 Ma in the Plio-Pleistocene [14].
The southern and northern slopes of the Himalayas differ dramatically in temperature and precipitation. Myricaria species occurring on the northern slope are generally xeric, same as those of the main plateau [73]. The western portions of the Himalayas and adjacent QTP are more arid than the eastern parts [27,73], see Figure 1c. Probably these differences are the cause of the persistent diversification between the eastern (A) and western Himalayas (B) for the Myricaria lineages (Figures 1 and 3). The eastern A contains the alpine species M. rosea, and M. wardii, and the western B includes M. elegans and M. prostrata. These four species are endemic to the Himalayas, and occupy two of the four Myricaria sections. In particular, the western Himalayas (B) is an important geographical node and dispersal center for Myricaria, with movement toward the Pamir-Alai, Hendukosh, Tianshan, and Kunlun-Altun mountains, etc. Noticeable, the Himalayas as a union was divided into A and B two times (see Figure 3), once at the time of generic origin and diversification, and the other at the diversification node of the sections Renantherae and Laxiflorae. Overall, the Himalayan areas AB, A, and B as the ancestors occurred at least seven times at the nodes in Figure 3. In detail, the Himalayan union AB as an ancestral area appeared at two nodes, with estimated ages of 20.25 Ma (genus crown age) and 10.07 Ma, the eastern Himalayas A at three nodes with ages of 19.84 Ma, 6.83 Ma (section Alpinae crown age), and 6.35 Ma (section Laxiflorae crown age), and the western Himalayas B occupied two, with ages of 8.83 Ma (section Parallelantherae crown age) and 5.81 Ma. All of these stressed the diversification and geographical heritage of Himalayan ancestry.
In some cases, the Himalayas are regarded as a corridor of plant species migration linking East Asia and Central Asia and/or the Mediterranean [16,35,[74][75][76][77]. The Himalayas have acted as a migration corridor for Sino-Japanese elements westward and Mediterranean elements eastward, such as Hippophae rhamnoides (Elaeagnaceae) [37], Pogonophace (Astragalus, Fabaceae) [72], and Phyllolobium [77] (Fabaceae). They also served as a migration route westwards to Central Asia from East Asia for Caragana (Fabaceae) [78], Begonia migration was through the Himalayas eastwards from Africa [40]. Our Myricaria analysis shows endemism, origin, and remarkable multiple diversifications within the Himalayas, but with weak or absent migration here. The Himalayas therefore seems to serve as the center of origin for Myricaria rather than as a corridor. This is rather like as these QTP endemic plant groups, such as Nannoglottis (Asteraceae) [79], Ligularia complex (Asteraceae) [80], Saussurea (Asteraceae) [81], Dolomiaea, Diplazoptilon and Xanthopappus (Asteraceae) [82], and Aconitum (Ranunculaceae) [83], all of them are hypothesized to have originated from native QTP during Miocene-Pliocene.

Diversification of Myricaria sections
Our estimated crown ages for the four sections of Myricaria were about 8.83,6.35 Ma (Figure 2), this is a remarkable molecular clock response to the rapid and major-range uplift of QTP at ca. The time of Myricaria origin and diversification related to Himalayan uplift in early Miocene, and the subsequent sectional diversifications at about 8.83,6.35 Ma, correspond respectively to the two labels of the QTP uplift, Himalayan motion, and rapid and major-range uplift [21][22][23][24][25][26][27][28][29]. These diversification times of the genus and sections are very similar to those of the genus Caragana [38], which also has a 16 Ma time of generic origin and diversification ages of sections at about 8 Ma. However, here we provide the Himalayan places of origin and diversification for Myricaria, while those of Caragana remain unknown. The evolution of Myricaria is also temporally similar to that of the Asian spiny frogs Paini [42]. Uplift of the Himalayas, is hypothesized to be most likely due to a cut off of the genetic exchange at ca. 19 Ma, resulting in the splitting of the subgenera Nanorana and Paa of Paini began in the Miocene near 10 Ma [42].

Migration along Asian mountains
For plant migration and dispersal, mountains generally act as a route or corridor [84], such as the Himalayan corridor mentioned above. Whereas the Himalayas for Myricaria are regarded as the center of origin, other distributions outside of the Himalayas and QTP can be understood as dispersals or migrations eastward, westward and northward along Asian mountains (see Figure 1). In fact, the results of vicariance and dispersal from biogeographical analysis (Figure 3) show that except for divergence of phylogenetically basal clades located in the eastern and western Himalayas (A and B), most remaining events, occurring in areas such as the Tianshan-Jungger-Turan (C), Altai-Sibiria (D), Mongolian Plateau (E), and Asia Minor-Caucasus-Europe (F) resulted from dispersal events. As evidenced from Myricaria (see Figure 1), dispersal and migration were possible from the Himalayas to Asia Minor-Caucasus-Europe as shown by M. germanica, and to the sand lands of the Mongolian Plateau by M. bracteata. These distributions along the Asian mountains are very similar to those of Hippophae rhamnoides (Elaeagnaceae) [37], a species with another interesting distribution in North Temperate Eurasia. Hippophae rhamnoides includes nine subspecies, and has been shown to have originated from the QTP, or more exactly, the eastern QTP-Hengduan Mountains, and then to have radiated and dispersed in different directions. Here, Myricaria originated in the Himalayas union, not in the eastern Himalayas only as H. rhamnoides. However, the northwestern Himalayas for H. rhamnoides and Myricaria played an important node role in connecting with Central Asia and Europe, and both dispersal route and direction is very similar.

Supporting Information
Appendix S1 Two new sections within Myricaria. (DOC)