Genetic Diversity and Evolutionary Relationships of Chinese Pepper Based on nrDNA Markers

Chinese pepper, referring to Zanthoxylum bungeanum Maxim. and Zanthoxylum armatum DC. species, is an important spice crop that has long attracted people’s interest due to its extensive application in Asian cuisine to improve taste. Numerous cultivars have been developed during the long history of domestication and cultivation. However, little to no information is available on the genetic diversity and evolutionary relationships of Chinese pepper cultivars and their historical diversification has not been clarified. Herein, we sequenced two nrDNA non-coding region markers, the external transcribed spacer (ETS) and the internal transcribed spacer 2 (ITS2), to assess genetic diversity and phylogenetic relationships among 39 cultivated and wild populations of Chinese pepper from eight provinces and to address the question of ancient demographic trends which were probably influenced by changing climate during evolutionary history. In total, 31 haplotypes were identified based on 101 polymorphism sites. Our results revealed relatively high level of genetic variation despite long-term cultivation of this crop. AMOVA revealed that genetic variation existed predominantly within provinces rather than among provinces. The genetic structure result based on haplotype network analysis largely reflected historical records, which suggested a Gansu origin for Chinese pepper and an ancient west-to-east spread of Chinese pepper circulating in China. We also provided evidence that changing Pleistocene climates had shaped the demographic trends of Chinese pepper. Taken together, our findings not only suggest that Chinese pepper is a dynamic genetic system that responds to evolutionary forces, but it also provides a fundamental genetic profile for the conservation and responsible exploitation of the extant germplasm of Chinese pepper and for improving the genetic basis for breeding the cultivars.


Introduction
Zanthoxylum, a genus of both economical and ecological importance in the Rutaceae family, consists of 250 species spreading in tropical and subtropical areas, including Asia, Africa, Oceania, and North America [1]. To date, at least 45 species and 13 varieties in Zanthoxylum have been found within China, representing one of the modern diversity centers of the genus [2]. At least two species within this genus, commonly referring to Z. armatum and Z. bungeanum (collectively referred to as "Chinese pepper" henceforth), are used in Chinese cuisines to improve the taste, mainly due to their numbing and tingling oral sensation created by a chemical component called alkylamides [3][4][5]. Chinese pepper has been demonstrated to possess a wide range of definite anti-tumor and pharmacological effects on, i.e., digestive system and nervous system. These salubrious effects appear to be most likely due to analgesic, stress-reducing, anti-inflammatory, and antioxidant effects [9,10]. In addition, these two species have ecological values and are planted as mountainous afforestation trees because of their ability to conserve soil and water [11]. Therefore, both Z. bungeanum and Z. armatum are excellent tree species for cultivation with economical, medicinal, and ecological values.
During the long history of domestication and cultivation of Chinese pepper, numerous excellent local cultivars have been developed. However, the genetic diversity and evolutionary relationship among these cultivars remain unexploited. There is still lots of confusion in the classification of cultivars. Traditionally, the identification of cultivars are mainly based on color, shape, fragrance, and mature period, which possibly led to lots of critical problems in classification. For example, the 'synonym' or 'homonym' phenomenon is prevalent, leading to inconvenience for scientific research and market. Our previous study has shown that some cultivars (i.e., Qinghuajiao) thought to belong to Z. schinifolium are in fact nested within Z. armatum, which indicated that earlier cultivar classifications need revision [12]. Therefore, in order to adequately preserve and utilize the germplasm resources of Chinese pepper, it is necessary to understand the genetic diversity, genetic structure and evolutionary relationship of cultivars using more effective and specific methods, such as DNA markers.
Recently, the genetic diversity of some Zanthoxylum species have been estimated based on random-amplified polymorphic DNA (RAPD) [13,14], intersimple sequence repeat (ISSR) [13,15], and amplified fragment length polymorphism (AFLP) markers [16]. However, these dominant markers neither discriminate between homozygotes and heterozygotes, nor provide accurate estimates of allele frequencies. Furthermore, these markers are generally unable to resolve more distant evolutionary relationships. Therefore, due to the high level of polymorphism, codominance nature, biparental inheritance, and reproducibility, a mass of SSR markers were developed for Zanthoxylum species [17][18][19]. DNA sequence variation has been are most suitable for inferring evolutionary histories, and SNPs are a rapidly developing source of sequence data from multiple loci. Yoshida et al. infered the demographic history of Z. ailanthoides Sieb. et. Zucc. [20], a species native to Japan, using nucleotide variation. This research revealed a relatively recent divergence between these two current populations and rapid formation of strong genetic structure among subpopulations. These results were not clarified in their previous study using microsatellite markers [21]. Lee et al. performed a comparative analysis on the complete chloroplast genome sequences and metabolic profiling between Korean and Chinese Zanthoxylum species, in which they developed seven validated markers for distinguishing Z. schinifolium Sieb. et Zucc., Z. piperitum (L.) DC., and Z. bungeanum [22]. Although the gradual advances in genetic diversity and phylogenetic relationships over time, many questions remain unresolved as to the origins and evolutionary dynamics of Chinese pepper.
As a result of high accuracy, rapid evolution, and independence of environmental variation, the nrDNA non-coding region markers have been chosen for genetic studies. The internal transcribed spacer (ITS) regions, especially ITS2, are thought to be effective phylogenetic markers for identifying medicinal samples [23,24]. Similarly, the external transcribed spacer (ETS) regions are also commonly used in crop breeding system [25]. Here, we took advantage of our extensive sampling of cultivated and wild Chinese pepper populations occurring in eight provinces of China, and ETS and ITS2 markers, to achieve the following objectives: (1) estimate the genetic diversity of Chinese pepper; (2) determine the evolutionary relationships of Chinese pepper cultivars; and (3) resolve the domestication origin and spread routes of Chinese pepper compared to the written records.

Plant Materials
Young leaf materials were collected from healthy trees grown in eight main producing provinces, including Shaanxi, Henan, Shandong, Shanxi, Sichuan, Yunnan, Gansu, and Guizhou in China, covering 37 cultivars and two wild populations. Fresh leaves were dried immediately using silica gel dessicants (HG/T2765.4-2005, Qingdao Haiyang Chemical Co., Ltd., China) in the field. Three to five leaves were placed in a valve bag (12 × 17 cm) with one-third silica gel dessicants. These samples included 84 accessions of Z. bungeanum and 21 accessions of Z. armatum (Table S1 and Figure 2A). Particularly, we focused our sampling on accessions that have been locally cultivated for a long time and are at least one kilometer apart from each other to avoid clone.

Data Analyses
We used BLAST to examine the concordance of each site and determine the homology between these sequences and other related taxa sequences taken from GeneBank. All nrDNA sequences of Chinese pepper were deposited in NCBI under accession numbers of MK922480-MK922510. We chose the same nrDNA non-coding region sequences (downloaded from NCBI) of Z. kauaense and Z. ovalifolium as outgroups with accession numbers of HG002518 (ITS2) and MG975154 (ETS) as well as MH016527 (ITS2) and MG975167 (ETS). Sequences of haplotypes and two outgroups were aligned and concatenated using MEGA7 [28].

Genetic Diversity
Haplotype diversity (H d ) and nucleotide diversity (P i ) were calculated in DnaSP 6 with the nucleotide sequence chosen as DNA, diploid, and autosome [29]. In order to examine the genetic differentiation, analyses of molecular variances (AMOVA) were performed in ARLEQUIN [30].

Genetic Structure and Evolutionary Analyses
The pattern of genetic structure and evolutionary relationships among haplotypes were assessed using three different methods. The evolutionary relationship of 31 haplotypes was inferred using the BEAST 1.8 [31]. In order to initialize BEAST, jModelTest and the Akaike information criterion were used to obtain the appropriate nucleotide substitution models, favoring a GTR+I+G model for each alignment. We tested whether a molecular clock could be fitted to our data by comparing the ML value with and without the molecular clock constraints implemented in MEGA7. The null hypothesis of equal evolutionary rate throughout the tree (strict molecular clock) was strongly rejected (with clock, lnL = −4890.756; without clock, lnL = −3276.201; P < 0.0001). Thus, an uncorrelated log-normal relaxed-clock model was chosen. Thus, a Yule process and an uncorrelated log-normal relaxed-clock model were treated as tree prior. We conducted five independent BEAST runs with Markov Chain Monte Carlo (MCMC) simulations of 100 million iterations with sampling every 1000 generations, following a burn-in of the initial 20% cycles. The five runs were then combined in LogCombiner. The ESS of each parameter was greater than 200. The maximum clade credibility tree (MCC) was generated with TreeAnnotator 2.4 with the initial 10% of trees discarded as burn-in.
To determine the genealogical relationships of haplotypes, a haplotype network was constructed using TCS1.21 [32]. The TCS program collapses sequences into haplotypes and calculates the frequencies of the haplotypes in the samples. The statistical parsimony was calculated for pairwise differences until the probability was greater than 95%. In this analysis, indels were regarded as one single character resulted from a mutation event (i.e., treated as a fifth state). Finally, a neighbor-net tree was built from the concatenated dataset of 31 haplotypes with SplitsTree4 [33]. The parameters of split transformation, distance transformation, character transformation, variance, and bootstrap replicates were set to equal angle, neighbornet, uncorrected P, ordinary least squares, and 1000, respectively. Other parameters followed default settings.

Inference of Demographic History
In order to explore the demographic history of both Chinese pepper species, we examined the frequency distributions of pairwise mismatch between sequences using Arlequin 3.5 [30] and DnaSP [29]. The observed mismatch distribution was compared with that obtained under models of spatial expansion, population expansion, and constant population size for evidence of model fit by calculating the sum of squared deviations (SSD) of the observed data relative to the model and Harpending's raggedness statistic (H Rag ). To assess the significance of demographic expansion, we also used Arlequin and DnaSP to carry out Tajima's D, Fu's F s and Ramos-Onsins & Rozas's R 2 neutrality tests, which are the most powerful tests for detecting population growth with small sampling size [34].
To analyze changes in relative population sizes over time, we reconstructed demographic histories independently in both species by using the Bayesian skyline plot [35] as implemented in BEAST 1.8 [31], with the divergence time between Z. armatum and Z. bungeanum informed by previous phylogenetic analysis [36] as the calibration point. We modeled these calibrations as an normal distribution with a sigma of 1.4 and an offset (hard bound constraint) that equaled the median age of the calibration. This calibration point also generated an estimated mutation rate for Chinese pepper as 8.221 × 10 −9 substitutions per site per year. According to our long term field investigation, Chinese pepper have a generation time (the time needed to generate the next generation) of four years. Meanwhile, to obtain a more resolved picture on diversification histories, we performed a lineage-through-time (LTT) analysis with BEAST 1.8 to explore the temporal variation in lineage richness as a function of diversification rate within each species of Chinese pepper. The BEAST operations followed the parameters in Section 2.3.2, except choosing the Coalescent Bayesion Skyline as the tree prior. The final result was visualized in Tracer version1.6 to assess convergence and stationarity of each chain to the posterior distribution and then to generate the Bayesian skyline plot and LTT plot.

Sequence Characteristics and Haplotype Distribution
We surveyed 105 accessions of Z. bungeanum and Z. armatum. The ETS region evolved at a slightly faster rate than ITS2 with a higher variation ratio of 15.05% (62/412), than that of ITS2, 10.70% (46/430). After combination, the length of aligned nrDNA sequences was 845 bp with 98 substitutions and 3 indels. Based on 101 polymorphism sites, 31 haplotypes were defined and referred to as H1-H31, of which 28 (H1-H28) were specific to Z. bungeanum and three (H29-H31) to Z. armatum. Of the 31 haplotypes, nearly half (14/31) were unique (haplotypes that represented by only one accession); only five haplotypes represented by more than one accession were found in multiple provinces, and the remaining 12 were always confined to a single region. For Z. bungeanum, H1-H27 were only assigned to cultivars, while H28 corresponded to wild accessions. The H7 haplotype, the most abundant, included 25 out of 84 Z. bungeanum accessions that from Shaanxi, Sichuan, Yunnan, and Guizhou provinces. H4 accounted for the second most predominant haplotype found mainly in Shaanxi, Sichuan, and Gansu provinces. Surprisingly, other haplotypes only occurred in a single province. A large number of regionally unique haplotypes indicated that most cultivars were probably developed directly from local wild populations while only the elite ones could have the opportunity to be introduced to elsewhere. For Z. armatum, haplotypes H29 and H30 were derived from cultivated accessions, while only H31 were found in wild accessions. However, all three haplotypes appeared to have multiple origin, i.e., both H29 and H30 were found in Shaanxi, Sichuan, and Guizhou provinces, and H31 was found in Yunnan, Sichuan, and Guizhou provinces (Table 1).  The sampling sites could be naturally divided into three regions, northwest, north, and southwest China, due to Qinling Mountain acting as an important geographic and climatic dividing line and thus a natural barrier to gene flow for most plant populations. Based on the distribution pattern of haplotypes in China observed above, we can conclude that haplotypes were relatively abundant in the northwest China (i.e., Shaanxi and Gansu provinces), while Chinese pepper accessions in north China had a relatively more simple haplotype composition. We also observed numerous regional specific haplotypes, including H15-H18 and H21-H28 that were specific to the northwest China, H19-H21 specific to the southwest China, and H2-H3, H8-H11, and H12-H14 unique to the north China ( Figure 2A). These findings suggested limited inter-regional gene flow for Z. bungeanum species. However, we observed no apparent spatial genetic structure for variation of Z. armatum species due to its low number of haplotypes.

Genetic Diversity and Genetic Differentiation
At species level, Z. bungeanum showed both higher haplotype and nucleotide diversity than Z. armatum (Table 1). At a province level, Shanxi and Henan harbored the highest haplotype diversity (H d = 1.0000) for Z. bungeanum, which is probably caused by the high amount of introduced cultivars from other areas, as told by local farmers. While the highest haplotype diversity (H d = 0.5000) for Z. armatum appeared in Guizhou and Sichuan provinces, which also showed the highest nucleotide diversity. However, it should be noted that the high haplotype diversity found in these provinces may be due to the limited sampling. It is not unlikely that with more sampling we would have observed more identical sequences occurring in these regions, which would certainly reduce their haplotype diversity. For Z. bungeanum, the highest level of nucleotide diversity was found in Gansu (P i = 0.0267). However, the Z. bungeanum accessions in Guizhou and Z. armatum accessions in Yunnan both had only one single haplotype, thus their haplotype and nucleotide diversity were both zero (Table 1).
AMOVA was performed to assess how intra-or inter-specific genetic variation was partitioned between species and among/within geographic localities. The AMOVA showed significant differentiation between the two species, with 51.61% of total variation was due to differentiation between species, 11.17% was ascribed to differentiation among provinces of a species, and 37.22% of the genetic variation was due to differentiation among individuals ( Table 2). The results of the AMOVA analysis confirmed that most genetic diversity of Chinese pepper was because of differentiation between species, with low differentiation among provinces within species. When the AMOVA was calculated separately for Z. bungeanum and Z. armatum, the majority of genetic variance of both species was explained within, rather than among, provinces.

Genetic Structure
The MCC tree from BEAST analysis showed that all haplotypes were separated into two highly divergent and strongly supported clades, Clade I and Clade II, corresponding to Z. bungeanum and Z. armatum, respectively ( Figure 3A). Z. bungeanum accessions could be further divided into two sub-clades, referred to as sub-clade I and sub-clade II. Sub-clade I consisted of 15 haplotypes (H3, H6,  H8, H9, H11, H13, H14, H17, H18, H22, H23, H24, H26, H27, H28), derived from Shaanxi, Shanxi, Gansu, Henan, and Shandong, namely the north and northwest China. Sub-clade II included 13 haplotypes (H1, H2, H4, H5, H7, H10, H12, H15, H16, H19, H20, H21, H25) that derived from Hennan, Shaanxi, Gansu, Sichuan, Yunnan, and Guizhou. Haplotypes of Z. armatum could be classified into two sub-clades corresponding to cultivated (Sub-clade III) and wild (Sub-clade V) accessions, respectively. The results indicated a relatively deep divergence between cultivated and wild populations after the long history of domestication. A neighbor-net tree was built from concatenated sequences of 31 haplotypes using SplitsTree4 ( Figure 3B). This approach identified the same divergence pattern of haplotypes as shown in the MCC tree, but they differed in their ability to resolve the detailed relationship of Sub-clade I and Sub-clade II haplotypes (Figure 3). The evolutionary relationships of Chinese pepper haplotypes were further investigated with the haplotype network, the topology of which was perfectly consistent with the MCC tree ( Figure 2B). Three Z. armatum haplotypes formed a single group despite the wild haplotype H31 having eight mutational steps to the cultivated haplotypes H29 and H30. Three main haplotype groups can be recognized for Z. bungeanum species which was clearly separated from Z. armatum with at least 16 mutational steps between them. In general, ancestral or founding haplotypes are usually located at the center node of the groups and are more frequently distributed, whereas new formed haplotypes during recent differentiation events were on the edge of the network and distributed in specific regions. Accordingly, however, no ancestral or founding haplotypes could be inferred in this network because the predominant haplotypes-i.e., H7, H4, and H28-did not locate in the centers of each group, which indicated that these prevalent haplotypes were newly developed and afterwards spread widely. If the hypothesized haplotype at which Z. bungeanum and Z. armatum split was considered as the ancestral haplotype (i.e., a medium vector mv1 in Figure 3B) for both species, we uncovered clues about the origins and directionality of ancient spread of Z. bungeanum in our study regions. H1 and H5 were the most closely relatives to mv1 with only one or two mutational steps. Interestingly, these two haplotypes were only observed in Gansu and Shaanxi, the inferred origin of cultivated Chinese pepper from history records and the home to 'Shujiao' and 'Qinjiao' (the name of Chinese pepper recorded in ancient literatures). Therefore, we can confidently assumed that H1 was the most recent ancestor of Z. bungeanum Group II while H5 was the most recent ancestor of H7 and Z. bungeanum Group I. In Group II, H21, H4, H25, H15, and H16 were all sympatric haplotypes of H1, indicating local diversification of this group. The direct connection between Yunnan haplotype H19 and H4 largely ruled out the probability of local evolution but provided a potential evidence for long-distance dispersal of this haplotype. Besides H7, the Henan haplotype H2 and Shangdong haplotype H6 also closely related to H5. To be more specific, Shandong haplotypes all derived directly from Gansu ones (i. e., H5→ H6, H24→ H8, H23→ H9), and three out of four haplotypes in Henan also have found their direct ancestor haplotypes in Gansu (H5→ H2, H24→ H3, H23→ H11) and the remaining one can be traced back to its root in Shaanxi (H7/H1→H10). Afterwards, the Henan haplotype H3 continued to march northeastward to Shanxi and generated haplotypes H13 and H14 in succession there. We also observed short-distance eastward transmission from Gansu to Shaanxi (H5→ H7, H21→ H4) and local evolving (H27→ H17, H27→ H26, H4→ H25→ H15→ H16, Gansu H21→ H4). However, some links between nodes in the network appeared to be indications that some haplotypes backflowed to Gansu or Shaanxi (i.e., H14→ H22, H14→ H23, H6→ H24, H13→ H18) ( Figure 3B). This eastward spread was consistent with the history records that the increasing trade of Chinese pepper has activated the introduction of 'Shujiao' seedlings to Shandong.

Demographic Trends
In order to examine the impact of past climate change (i.e., glacial-interglacial cycles) and anthropogenic activities (i.e., domestication and transmission) on demographic trends of Chinese pepper, we performed several population dynamics inferences. The mismatch distribution for Z. bungeanum resembled a bell curve with a mean of 18.566 mismatches (Table 3) and a tail having a high frequency of pairs with only a few mismatches ( Figure 4A). This observed mismatch distribution did not differ significantly from the modeled distribution for either a stable population size or a spatial expansion model (H Rag , P = 0.948 and P = 0.247, respectively) but did differ significantly from the expectation of a sudden population expansion (P = 0.002, Table 3). However, we observed a pattern of better fit of the observed data to the spatial expansion model by a SSD nearly two times higher for the constant population size model (SSD = 0.025) compared with the spatial expansion model (SSD = 0.014). This finding indicated that the initial range of Chinese pepper population largely increased over time and over space and became generally subdivided in the sense that individuals would tend to mate with geographically close individuals rather than remote individuals. For species Z. armatum, all demographic models were rejected by mismatch distribution (SSD P < 0.05 for all three models), indicative of a sudden population contraction hypothesis ( Table 3). The three neutrality test statistics-Tajima's D, Fu's Fs, and Ramos-Onsins and Rozas's R 2 -all yield non-significant positive values for both species, largely indicative of neutral selection or population evolving as per mutation-drift equilibrium. However, a more powerful approach for population history inference, the Bayesian skyline coalescent model, suggests that both species has undergone a sudden recent decline in effective population size, despite that the lineage through time (LTT) plots indicated that both Chinese pepper species did not exhibit a constant diversification rate: an acceleration in diversification of the Chinese pepper lineages occurred in the Pleistocene (ESS > 200) (Figure 4).

Genetic Diversity and Genetic Differentiation
The genetic diversity within a species/population is the outcome of its selective and demographic past and, in turn, it determines its evolutionary potential and thus long-term fitness and adaptation. Genetic differences between populations are presumed to be a consequence of different genomic DNA sequences following mutation and subsequent selection, drift, and other factors that alter allelic frequencies. Therefore, investigating genetic diversity could not only help in the understanding the evolution history of species but also managing genetic resources of important crops [37]. Chinese pepper is a specialty spice in China, where it has been bred for more than 2000 years, thus a large number of elite cultivars have been widely named, released, and dispersed. However, the evolutionary relationships and genetic diversity of Chinese pepper are still unresolved. Here, we have examined the genetic diversity and evolutionary relationships of Chinese pepper using the most comprehensive Chinese pepper sample set compiled to date to our knowledge, encompassing 37 local cultivars collected from the main growth area and two wild populations.
In general, the nucleotide diversity and haplotype diversity in most crops should have decreased dramatically for cultivated populations compared with their wild ancestors. This is because the strong selections for specific traits, combined with a series of population bottlenecks, severe or mild, during domestication have greatly alter the genetic structure of populations and the underlying genetic architecture of phenotypic traits. Often only 20% of the total diversity contained within the wild ancestor can be maintained through domestication [38]. However, we observed a increasing trend instead of decline in genetic diversity of cultivated population in both Z. armatum and Z. bungeanum species. Nevertheless, we note that this does not preclude the possibility of opposite scenario because the limited sampling of wild populations has probably led to underestimate of their genetic diversity, and more data will need to be collected to address this issue.
Subsequently, we compared the genetic diversity difference between these two species. As shown, Z. armatum exhibited lower genetic and haplotype diversity. This result conflicts with the previous genetic analyses based on three cpDNA markers [12]. This is partially attributable to the difference in evolution history between cpDNA and nrDNA sequences. The set of nrDNA regions that an offspring inherits from both parents through a combination of the genetic material of each, certainly have relatively rapid rates of nucleotide substitution when compared with the maternally inherited cpDNA [39].
We observed that Gansu had the highest P i for Z. bungeanum, which confirmed our assumption that Gansu was once one of the origin or genetic diversity centers of Chinese pepper. Interestingly, this finding, combined with the eastward migration of some haplotypes in Gansu and Shaanxi (directly to Shandong), have largely reflected the history records, such as in the book Qimin yaoshu. The evolutionary network of haplotypes was an efficient approach that can be used to recover the directions of haplotypes transmissions at the local scale and spread at the larger scale, and particularly, to identify the central haplotype that can be considered as the ancestral or founding haplotype. The evolutionary network of haplotypes in this study indicated that Gansu deservedly harbored the same sequence types as H14, H6, and H13 that were not contained in our dataset. More sampling in this region will give us a good chance to compare sequences among these accessions. Thus, we suppose that Gansu haplotypes H22, H23, and H24 and Shaanxi haplotype H18 arose from nucleotide mutations of local counterparts of H14, H6, and H13, respectively.
The most frequent haplotypes, designated H7 and H4, were each represented by 25 and 13 samples that emerged in several provinces. So where did H7 and H4 originated? H7 derived from Gansu-specific H5 by six mutational steps and occurred in Shaanxi, Sichuan, Guizhou, and Yunnan. Given the geographic proximity, H7 most probably appeared first in Shaanxi and then was introduced southward by farmers due to its excellence in quality and yields. H4 occurred in Gansu, Shaanxi, and Sichuan and was derived directly from H21, that originated in Gansu and Sichuan. Based on the tenet of nearby domestication (domesticating from wild populations nearby), the most probable scenario was that H4 evolved from local H21, despite H21 sequence type was not observed in Shaanxi, but we were sure of the existence of H21 in Shaanxi because its direct ancestor H1 was found there.
As shown by AMOVA, Chinese pepper harbored substantial within province and intraspecific variation compared with interspecific variation, which may derive from several evolution processes, including high gene flow among provinces, microvicariance and subsequent population mixing, and so on. Despite this, we did observe geographic isolation effects that generated both genetic differentiation and phenotypic diversity around, i.e., Qinling Mountains. The two most elite cultivars in Shaanxi, FXDHP and HCDHP, each represented by H4 haplotype and H7 haplotype, were located, respectively, north and south of Qinling Mountains. Furthermore, FXDHP and HCDHP showed great differences in terms of morphological characters and quality. The leaves of HCDHP are narrow but thick with light green color, whereas FXDHP leaves are wider, more flat and dark green. The pericarp of HCDHP is thicker with slightly yellow color, having fewer, but bigger, oil glands. The fruits of HCDHP imparted a persistent numbing sensation and a light aroma flavor, but with a slightly bitter aftertaste. However, the FXDHP has smaller fruit size, with well-developed oil glands on the pericarp, and engender a stronger numbing and tingling feelings in the tongue without any bitter aftertaste. These findings may have suggested that both natural processes such as isolation, local adaptation, and drift and the agriculture practices of local farmers such as selection, seed exchange, have largely shaped the genetic structure and morphological diversity of Chinese pepper.

Demographic Trends Were Largely Shaped by Pleistocene Climate Oscillations
A genetic signature of population contraction in both Chinese pepper species was reflected in the coalescent Bayesian skyline plot, which also allow for estimates of the timescale of ancient contraction event. Provided a generation time of four years and that the split between Z. bungeanum and Z. armatum occurred 15 million years ago (Mya), we can conclude that both species showed signatures of population decline during late Pleistocene (i.e., 1 Mya to present). Interestingly, this timescale coincided with that of the rapid diversification in both species. Thus, they were likely accordant with time periods that correspond with known events of environmental upheaval. During this time, subtropical China experienced three glacial periods (i.e., Naynayxungla Glaciation, Penultimate Glaciation, and the last glacial maximum) and two interglacial periods, suggestive of crucial phases of evolution of the East Asian Pleistocene glaciation [40]. Based on our results of demographic change and LTT plots, we explored the possibility that drastic environmental changes during this period has induced the fragmentation and variation of ancestral suitable habitat and, thus radiation of the Chinese pepper species. We proposed that, during the late Pleistocene climate fluctuations, the ancestor populations of Chinese pepper species might have been forced into ongoing cycles of retreating into, and the re-expansion from, refugia. This was confirmed by the finding of recent spatial expansion of Chinese pepper in the sequence mismatch distributions, despite that the exact timescale of expansion was difficult to obtain due to a lacking of known mutation rate for Chinese pepper non-coding regions. However, we could obtain a rough estimate of the timing of spatial expansion event based on a mutation rate of 8.221 × 10 −9 substitutions per site per year derived from BEAST analysis, which fell within a common believed range in angiosperm plants, 5.0-30.0 × 10 −9 substitutions per site per year [41]. This mutation rate generated an estimate of time, t = 1 million years, since expansion, according to the formula τ = 2µt [42] (µ is the mutation rate per nucleotide per generation derived from BEAST analysis), just overlapping with the Pleistocene cycle of glaciation. Thus, we hypothesized that recurrent episodes of harsh climate conditions have resulted in the isolation of small populations by mountains over generations, which shrank the populations but expedited diversification and the fixation of morphological traits. This process was probably reinforced during the recent thousands of years of local domestication by farmers. The two species showed similar trends of effective population size drop but the rate of diversification increased during the glacial periods, indicating strong robustness and reliability of our inference.
It is worth noting that structured populations generate signals of population size change [43,44]. Therefore, it is always difficult to determine whether demographic events such as expansions or contractions (bottlenecks) inferred from genetic data are real or due to the fact that populations are structured in nature. If two genomic haplotypes sampled derive from the same subpopulation when sampling, the inferred demographic trend depicts a decrease of the effective population size similar to a bottleneck signal. Conversely, if each of the two haplotypes sampled derive from different subpopulations, the resulting inferred demographic trend is the population expansion. In our study, the sampling covered the whole distribution region of Chinese pepper in China, which likely could reduce the bias of demographic inference derived from structure population.
What about the robustness of our inference based on a small set of genetic markers? Understanding the relative role of neutral and selective forces in shaping the overall genetic diversity has been one primary aim of population genetics. This is often achieved by investigating the distribution patterns of genome across multiple samples. It is well known that the genotyped dataset of variation sites may suffer from an ascertainment bias due to the procedure used to detect variations. The sequencing of a set of gene regions is a common SNP genotyping methods recently used. However, such SNPs tend to be of higher frequency than random SNPs scattered across the whole genome and may not be geographically representative. Additionally, the genes in the genome are not evolving at the same rate, the genetic regions under selection pressure in general preserve lower degree of genetic variation than fast-evolving non-coding genome regions. Therefore, such ascertainment biases probably result in biased frequency spectrum toward common alleles, and various factors, such as population subdivision, recombination, and measures of variability may be biased. For example, our sequencing of a small set of genetic regions lacks the potential of discovering rare alleles, which in turn may probably lead to an overestimate of Tajima's D.
In the context of a small set of loci, one important aspect of designing a successful study is to determine the appropriate number of samples to collect and analyze. The sample size must be large enough to ensure that the results are shaped by actual ecological patterns rather than by chance. For analysis of intra-species or population-level differentiation, fast-evolving genes or DNA regions are most suited. For our analysis, our sampling covered the main distribution regions of Chinese pepper and are enough to represent the genetic resources of Chinese pepper. Our chosen nrDNA markers are neutrally evolving and have enough variation sites that can be used to detect a significant difference between populations.

Conclusions
Chinese pepper includes at least two spice species of ecological and economical importance in China which has been excessively exploited since 2000 years ago. Many ecological and geographical variants of Chinese pepper have been formed in the process of long-term effects of the geographic environment and human domestication. In order to protect and manage the genetic resources of these important species efficiently, greater knowledge of the genetic diversity and population structure is needed. In the current study, we assessed the genetic diversity and population structure of eight provinces covering 37 cultivars and two wild populations representing most of the Chinese pepper range, by using two nuclear sequences. First, we found that the level of genetic diversity in the cultivated populations was still high, and the level of genetic differentiation among provinces was very low. Second, our results provided evidence for an ancient west-to-east spread of Chinese pepper, which is of particular interest in light of both ancient records and our evolutionary network which suggested a Gansu origin for Chinese pepper circulating in China. Additionally, we found evidence for the impact of changing Pleistocene climates on the demographic trends of Chinese pepper. Together, our findings will not only inform efforts for the conservation and management of the genetic resources of Chinese pepper, but also provide guidance for future studies of population genetics and breeding programs.
Author Contributions: Conceptualization, A.W. and S.F.; Methodology and software, Z.L.; Investigation and resources, S.F., L.T., and X.W.; Data curation, S.F.; Writing-original draft preparation, J.N. and S.F.; Writing-review and editing, Z.L.; Project administration, A.W.; Funding acquisition, A.W. All authors have read and agreed to the published version of the manuscript.
Funding: This research was financially supported by the National Key R&D Program of China (2018YFD1000605).