Phylogeography of two closely related species of Allium endemic to East Asia: Population evolution in response to climate oscillations

Abstract This study investigated the effects of climate oscillations on the evolution of two closely related Allium species, A. neriniflorum and A. tubiflorum. We sequenced three cp DNA (cpDNA) fragments (rps16, rpl32‐trnL, and trnD‐trnT, together approximately 2,500 bp in length) of two closely related Allium species, with samples from 367 individuals in 47 populations distributed across the total range of these species. The interspecific and intraspecific divergence times of the two species were in the Quaternary glaciation. The population divergence was high for the cpDNA variation, suggesting a significant phylogeographic structure (NST = 0.844, GST = 0.798, p < 0.05). Remarkable ecological differentiation was also revealed by Niche models and statistical analyses. Our results suggest the speciation event of the two species was triggered by violent climatic changes during the Quaternary glaciation.


| INTRODUC TI ON
Climatic oscillations and related changes in vegetation can be recognized in the period from 3.6 to 0.8 Ma (Naidina & Richards, 2016).
Multiple studies have showed that the climate oscillations especially in the glaciation cycles of the Quaternary may have facilitated speciation, diversification, hybridization, and changes in the distribution of the global vegetation (Liu et al., 2018;Shahzad, Jia, Chen, Zeb, & Li, 2017;Yamamoto, Ohtani, Kurata, & Setoguchi, 2017).
Although East Asia was profoundly affected by global climatic oscillations, even in the Quaternary ice ages, no massive ice sheet development occurred in this region (Hewitt, 2000). Therefore, East Asia was the most important refuge for many species during these climate fluctuations (Qiu, Fu, & Comes, 2011). Fossil-based biome reconstructions of East Asian species predict that, during the last glacial maximum (LGM), the temperate forests that currently dominate East Asia retreated southward to 25°N to 30°N (Harrison, Yu, Takahara, & Prentice, 2001). Nevertheless, recent phylogeographic studies have highlighted in situ survival of hardy alpine herbs and forest trees during glacial periods in northern China. Many refugia scenarios have been revealed for these species, such as single refugium (Liu & Harada, 2014), multiple refugia (Tian et al., 2009), refugia within refugia (Wang, Gao, Kang, Lowe, & Huang, 2009), cryptic refugia, or microrefugia (Stewart, Lister, Barnes, & Dalen, 2010).
Orographic influences have also had a crucial role in the interintraspecific divergence and population demography of speciation.
Recently, many phylogeographic studies have assumed a causal link between geological processes (orogenesis) and biological responses (diversification) (Chen et al., 2012;Jiang, Zhang, Zhang, & Sanderson, 2014;Wang & Yan, 2014;Wang, Zhang, & Yin, 2016). The results According to flora of china, we went to collect samples of three species in Sect. Caloscordum in the field. However, during our field sampling, we just found A. neriniflorum and A. tubiflorum and did not find A. inutile. Previous studies investigating A. neriniflorum and A. tubiflorum have shown that the two species are sister groups and have a close relationship (Li & Xu, 1996;Lu, Yang, Lu, Zhou, & He, 2016). The leaf micromorphology characters of A. neriniflorum and A. tubiflorum are different, which can serve as a classification basis (Lu et al., 2016).
In this study, we sampled populations across the distribution range of A. neriniflorum and A. tubiflorum for phylogeographic analyses and used three sets of cpDNA fragment sequences to construct phylogeographic patterns for these two species. We aimed to answer the following three questions: (a) Are A. neriniflorum and A. tubiflorum significantly different from molecular levels? Is the genetic diversity of A. neriniflorum and A. tubiflorum analogous? (b) Where were the refugia of both species during the climatic oscillations? Did they survive in situ or retreat to somewhere below to 30°N during climatic oscillations? (c) How did the climatic oscillations or other geological factors effect on the population evolutionary history of these two species?

| Population sampling
According to flora of china, we went to collect samples of three species in Sect. Caloscordum in the field. However, during our field sampling, we did not find A. inutile in Anhui, China. As for A. neriniflorum and A. tubiflorum, we found that their morphological characters are very similar. According to our statistical analysis, the main difference in the shape of these two species based on quantitative characters is that the plant height and pedicel length, A. neriniflorum are usually taller in plant height and longer in pedicel length than A. tubiflorum (see Figure 1). We also found that the ecological conditions of A. neriniflorum and A. tubiflorum differ. A. neriniflorum is easily found in damp places and meadows of low altitude. In contrast, A. tubiflorum mainly grows in the   (Table 1, Figure 2). Voucher specimens (Table 1) were also collected for all sampled populations. The specimens were deposited in the archives of the Herbarium of Sichuan University (SZ).

| Network representation
The genealogical relationships among all the cpDNA haplotypes were examined via a median-joining (MJ) haplotype network, which was constructed using NETWORK ver. 5.0 (Bandelt, Forster, & Rohl, 1999). In addition, to complement the results of NETWORK 5.0, we used TCS 1.21 (Clement, Posada, & Crandall, 2000) to construct haplotype relationships (Templeton, 1995). Our analysis indicated that both site mutations and indels were equally likely to evolve and that each indel originated independently of other indels.
F I G U R E 2 Map of the sampling locations and the geographic locality of the cpDNA haplotypes detected for Allium tubiflorum and Allium neriniflorum.
Black circles indicate the distribution of A. neriniflorum populations; blue circles represent the A. tubiflorum populations.
In addition, the white circles represent probable hybrid populations. The cpDNA haplotypes (H1-H25) and their frequencies in each population are indicated by the colorful pie charts. The population codes and information are shown in Table 1

| Genetic structure of populations and geographical structure analyses
The number of haplotypes, haplotype diversity (H d ) and nucleotide diversity (π), were estimated with DNASP version 6.0 (Rozas et al., 2017). The average gene diversity within populations (H S ), total gene diversity (H T ), and two-population differentiation parameters (G ST and N ST ) at the species level were statistically compared using 1,000 permutations in the program Permut CpSSR 2.0 (Pons & Petit, 1996). Analyses of molecular variance (AMOVA) was applied in the program Arlequin ver 3.5 (Excoffier & Lischer, 2010) to assess the genetic differentiation within and between populations, while testing for the significance of partitions with 1,000 random permutations. The gene flow N m among populations was obtained from the following formula: F ST = 1/(1 + 2N m ) (Wright, 1950).

| Demographic history analyses
To detect whether the population groups (all populations, A. tubiflorum and A. neriniflorum) of A. tubiflorum and its congener A. neriniflorum experienced a recent population expansion, a mismatch distribution analysis was carried out using the DNASP program (Rozas et al., 2017). A multimodal distribution of differences between the haplotypes in a given population suggests the population did not experience an expansion, whereas a unimodal distribution indicates the occurrence of an expansion. Additionally, Tajima's D (Tajima, 1989) and Fu & Li's D* were calculated (Fu, 1997) with DNASP.

| Ecological niche modeling
In addition to the sample sites from this study, 57 collection records were obtained from the Chinese Virtual Herbarium (http://www. In Maxent, due to the small sample size of our data (15-79 localities), only the linear + quadratic + hinge functional forms were used (Merow et al., 2013). Here, models were built using linear + quadratic + hinge features, random training data were set to 25%, the create response curve and do jackknife options were selected, the logistic output format was logistic, and other parameters were set by TA B L E 2 Results of analysis of molecular variance (AMOVA) for the 46 populations of Allium tubiflorum and Allium neriniflorum based on cpDNA data (rps16, rpl32-trnL and trnD-trnT) Wilks' λ*** 0.103 (χ 2 = 115.885) *, ** and ***Indicate the value is significant at p < 0.05, p < 0.01, and p < 0.001, respectively, and # indicates the value is not significant (p > 0.05) based on MANOVA.
default. The accuracy of the model performance was tested by the AUC of the receiver operating characteristic (ROC) (Fawcett, 2006); values above 0.7 indicate that an ecological niche modeling (ENM) is useful (Araújo, Pearson, Thuiller, & Erhard, 2005).

| Test for niche overlap
To visualize the geographical niche zone, the overlap where hybrids between A. tubiflorum and A. neriniflorum may occur, we calculated the overlapping area based on the binarized layer using ArcGIS 10.2 (ESRI, Redlands, CA, USA). We used a binary threshold from DIVA-GIS of zero, which excluded the zero-suitable habitats, and the minimum training presence cumulative threshold was from Maxent (Zhao, Gugger, Xia, & Li, 2016).

| Environmental variation between Allium tubiflorum and Allium neriniflorum
Three statistical approaches were used to investigate the ecologi- Finally, Wilks' k was used to test the null hypothesis that A. tubiflorum and A. neriniflorum had an equal mean across climate variables via discriminant function analysis (DFA) using SPSS. The significance of Wilks' k was tested by chi-square tests. If the p-value of chi-square is <0.05, a significant ecological difference between the two species can be concluded.

| Haplotype distribution and phylogenetic analyses
The aligned rps16, rpl32-trnL and trnD-trnT sequences of Allium were 819, 888 and 777 bp in length, respectively. A 2,500-bp data set with 36 nucleotide substitutions and 11 indels was obtained when all the three fragments were combined (Supporting Information   Table S1). All the sequences were deposited in the GenBank database, with accession numbers MG709262-MG709289 (rpl32-trnL), MG709290-MG709317 (rps16) and MG709318-MG709345 (trnD-trnT).  . These all suggest that they represent ancestor haplotypes (Ignazio & Kohn, 2001). However, haplotype H14 is very special and become another clade, which means that populations containing H14 maybe a hybrid species. The related research is in progress. AMOVA analysis revealed that, in all the sampled populations, the variation attributed to differences between species was 76.22%, which is very significant, while that variation among and within population differences were 15.85% and 7.93%, respectively ( Table 2). The F ST value obtained by AMOVA also indicated high differentiation within the populations (F ST = 0.92656, p < 0.001). Approximately 76% of the variation was attributed to the species differentiation, and among-population variation accounted for just one-quarter (16%) of the total variation (Table 3), which indicates that A. tubiflorum and A. neriniflorum are two different species. For A. tubiflorum, the among-population differentiation was 87%, higher than that of its congener A. neriniflorum (55%). This result showed that the gene flow among populations of A. tubiflorum was lower than that among populations of A. neriniflorum. Overall, these results strongly indicate that the haplotypes are geographically structured across the species distribution range.
Cp DNA data shows that A. tubiflorum had higher levels of genetic diversity based on cpDNA (H d = 0.826) than A. neriniflorum (H d = 0.734) (Table 1)

| Demography
To test whether population groups of A. neriniflorum and A. tubiflorum underwent a recent demographic population expansion, we plotted the mismatch distribution of each group based on the observed number of differences between pairs of haplotypes.
The mismatch distributions of the two species were clearly multimodal curves (Figure 4)

| Ecological niche modeling
The AUC value for the current potential distribution of A. neriniflorum and A. tubiflorum both were relatively high (AUC > 0.9), demonstrating a reliable predictive model performance. Model exhibited distinctive ecological niches between A. neriniflorum and A. tubiflorum (Figure 6a-c). PCA biplots displayed distinct habitats between A. neriniflorum and A. tubiflorum (Figure 6d). The bootstrapping MANOVA estimates indicated that most climate variables contributed significantly to niche divergence, except for the mean temperature of the warmest quarter and the precipitation of the wettest month (Table 3). Wilks' k indicated a significant climatic difference between the two species habitats (χ 2 = 115.885, p < 0.001) (Table 3).

| Genetic characteristics of populations of two species
Our cpDNA data support the validity of both species, A. neriniflorum and A. tubiflorum, although both species show close morphological similarities (Lu et al., 2016). Therefore, cp DNA can be used as a marker to distinguish them in plant taxonomy.
We also found some differences that A. tubiflorum had higher levels of genetic diversity and the lower gene flow compared to These less blocked ecological environments resulted in frequent gene flow among populations. Therefore, we can inferred that the mountain taxa showed a high biodiversity, which agrees with the "mountain-geo-biodiversity hypothesis". The third reason may be the differently configured features of species itself. We found that A. tubiflorum has 0.8-7 cm pedicels and 3-6 ovules per locule, while A. neriniflorum has 4.5-11 cm pedicels and 5-8 ovules per locule. From the above data, the shorter pedicels and lower number of ovules per locule of A. tubiflorum cause their seeds to spread a shorter distance than those of A. neriniflorum. The population genetic differences of monolepsis cpDNA reflect the pathways and distance of species seeds, as was also shown by the study of Gao and Ge (1999)

| Multiple refugia
Phylogeographic studies on cool temperate deciduous species, such as Ostryopsis davidiana (Tian et al., 2009), Mongolian oak (Zeng, Wang, Liao, Wang, & Zhang, 2015), and Juglans mandshurica (Bai, Liao, & Zhang, 2010), suggested the in situ survival of species and multiple refugia in the northern regions of modern distribution of temperate forests. As herbaceous plants sensitive to climate, did A. neriniflorum and A. tubiflorum survive in a similar way as the abovementioned species? The cpDNA phylogenetic relationship and haplotype network analyses of the two species indicated that they belong to two different lineages restricted to the northeastern and central areas of China (Figure 3). The genetic differentiation between these two areas was high. Moreover, mismatch distribution confirmed that populations from the two regions had not undergone recent expansions. These data are indicative of fragmentation events and the existence of at least two glacial refugia in the current range of the Alliums species (Aoki, Matsumura, Hattori, & Murakami, 2006;Chen et al., 2012).
Moreover, the network analysis shows that haplotypes H1 and H18 are ancestor haplotypes. If ancestral and derived haplotypes do not overlap and are located in different regions, then ancestral haplotypes should be found close to refugia, while derived haplotypes are more likely to occur at the leading edge of the range expansion (Hewitt, 2000). In A. tubiflorum, haplotypes H17 and H15 are distributed in the northeastern and southwestern leading edges, respectively, whereas the ancestral haplotype H18 is restricted to the six populations (P10, P12, P20, p25, p28, p30) in the edge of the Taihang Mountains and Qinling Mountains (Figure 2), the possible range of the refugium for A. tubiflorum during the glaciation.
This coincides with the sanctuary of temperate walnut tree and temperate deciduous shrub species (Ostryopsis davidiana Decne., Betulaceae) (Tian, 2007). For A. neriniflorum, our study found the ancestral haplotype H1 had a disjunctive distribution in the northernmost edge (population P14) and central region (population P49, P2) of the range (Figure 2). Moreover, the disjunctive distance was nearly 500 km. The proposal is that, in the past, haplotype H1 was widely dispersed in the range of the examined populations and later was replaced by other haplotypes in the middle areas of the studied range owing to the continuous contraction and colonization during the Quaternary glacial and interglacial periods. Under this proposal, population P14 in Lamadianzi (HeiLongjiang), a very damp grass lake with relatively steady climatic conditions, may represent a cryptic northern refugia during the glacial period (Li, Shao, Lu, Zhang, & Qiu, 2012). Additionally, the ENM analyses indicated that these regions were suitable distribution ranges for the two species, which also supported the possibility of these refugia.

| Evolution of two species
To investigate what paleoclimatic events triggered the divergence of A. neriniflorum and A. tubiflorum, the divergence time was estimated based on genetic distances ( Figure 5). Our results show a population divergence of A. neriniflorum and A. tubiflorum at approximately 2.2094 Ma, which is coincident with the onset of the Quaternary glaciation (approximately 3-2 Ma). In this period, Taihang Mountains nearly formed a S-N-striking uplifted mountain range within the central North China block. However, the divergence time is different from that in Hauenschild et al. (2017), which may result from the different analysis methods. Our analyses is as solid as theirs, and that even if our results would be slightly older, our interpretation remains similar and plausible.
Our ENM analysis suggests that A. neriniflorum and A. tubiflorum each occupy a distinct climate niche ( Figure 6). Niche differentiation can provide the preconditions for the adaptive divergence of fragmented populations and subsequent speciation (Abbott & Brennan, 2014;Hewitt, 1996;Zhao et al., 2016). Speciation can occur largely as an outcome of ecological differentiation maintaining morphological differences, even after substantial gene influx between related species (Mckinnon et al., 2004;Nagy, 1997;Zhao et al., 2016).
We inferred that the populations of ancestral taxon about A. tubiflorum and A. neriniflorum were once widely distributed in north and central China in the warm period of the Pliocene. Then, a cooling period occurred from the Pliocene to the Pleistocene. Especially in the Quaternary glaciation, the climate became colder and drier, and there was an obvious alternation between glacial-interglacial cycles (Crippa et al., 2016;Tian, Wei, Cai, Wang, & Li, 2018). In glacial epoch, the ancestral populations began to contract and were limited in some frag- Almost simultaneously, the intraspecific divergence of A. tubiflorum F I G U R E 6 Potential habitats and ecological divergence of Allium neriniflorum and Allium tubiflorum estimated by Maxent (a-c) and principal components analysis (PCA) (d). Black points show occurrences of A. neriniflorum and A. tubiflorum. The area under the curve (AUC) >0.90 indicates the model performed well. The ecological separation between A. neriniflorum (yellow) and A. tubiflorum (gray) are shown; violet indicates the possible hybrid zone (c). The variations represented by PC1 and PC2 are 58.51% and 31.68%, respectively; A. neriniflorum (red) and A. tubiflorum (violet) are separated by the PCs occurred at approximately 1.168 Ma. These data show that the simultaneous intraspecific differentiation of A. tubiflorum and A. neriniflorum resulted from common climate events in the Quaternary (Yang, Dong, & Lei, 2009). Six great glaciations occurred during the Quaternary glaciation in China (Cui et al., 2011). It was previously unclear which glaciations had the greatest impact on the current vegetation of A. tubiflorum and A. neriniflorum. The intraspecific divergence time of the two species coincides with the time of the first glaciation (approximately 1.2 Ma), in the late stage of the Early Pleistocene (Cui et al., 2011). Therefore, it is most likely that the divergence of the main lineages of A. tubiflorum and A. neriniflorum can be attributed to the climate oscillation of the largest glaciation in the Quaternary (Zhang, Fengquan, & Jianmin, 2000;Zheng, Xu, & Shen, 2002). In a word, the evolution of A. neriniflorum and A. tubiflorum mainly resulted from the violent climate changes that occurred in the Quaternary glaciation.

ACK N OWLED G M ENTS
We are grateful to the editor and two anonymous referees for their constructive comments on the manuscript. We also thank the edi-

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
JTY and XJH designed the experiment. JTY collected the field data. JTY, SDZ, and DQH analyzed the data. JTY drafted the text. All authors revised it critically and significantly improved the text.