1 Introduction

Artemisia halodendron Turcz. ex Bess. (Asteraceae, Anthemideae, Subgen. Dracunculus) is one of the most common semishrub species in the Horqin Sandy Land in Northeast China. It is important for vegetation rehabilitation in Horqin Sandy Land because of its high ecological value, including that (1) it is the key species of the plant communities and landscapes studies in Horqin Sandy Land (Li, 1991); (2) it plays a key role in the vegetation-restoration process due to its high drought tolerance, anti-wind erosion properties, and sand-burial resistance (Dong et al., 2000 ; Li et al., 2002 ; Zhao et al., 2006 ). Previous studies on A. halodendron focused on aspects of population-distribution patterns (Chao et al., 1999 ; Cao et al., 2008 ), biomass allocation (Li et al., 2005 ), breeding distribution (Li et al., 2005 ), morphological characteristics and physiological adaptations (Zhou et al., 1999 ), root longevity (Huang et al., 2009 ), genetic diversity (Huang et al., 2011 , 2014), and establishment (Li et al., 2002 ) in Horqin Sandy Land. However, systematic comparison of genetic relationships of A. halodendron among populations from different hydrothermal regions has not yet been reported.

Horqin Sandy Land is located in the agropastoral transitional zone between the Inner Mongolian Plateau and the Northeast Plains (42°41′N–45°45′N, 118°35′E–123°30′E) and is one of the four largest sandy areas in northern China; it covers an area of approximately 139,300 km², of which up to 71,884 km² is desertified sandy land (Wang, 2003; Zhao et al., 2003 ). Landscape in this area is characterized by sand dunes that alternate with gently undulating lowland areas (Li et al., 2005 ). This area belongs to the continental semi-arid monsoon climate and is in the temperate zone, with a mean annual temperature (AMT) of 3–7 °C and mean annual rainfall (AP) of 350–500 mm (Zhao et al., 2003 ). Over recent decades, this region has undergone severe desertification (Li et al., 2000 , 2004) and has displayed the northern-moving phenomenon of the interlocked agropasturing area of North China in the most recent hundred years (Zhao et al., 2000 , 2002).

In the present study, we used chloroplast DNA (cpDNA) trnL–F to examine the genetic diversity ofA. halodendron. We specifically aimed to address the following questions: (1) What is the level of nucleotide diversity in A. halodendron from different hydrothermal regions? (2) How are the identified haplotypes distributed within and among populations? In particular, is there subdivision in the different hydrothermal-level populations (according to hydrothermal synthesis index)? We attempted to interpret the results to provide baseline genetic information pertinent to the restoration and management of degraded ecosystems in arid and semi-arid areas.

2 Materials and methods 2.1 Sampling

A total of 243 individuals were sampled from 10 natural A. halodendron populations (Tables 1 and 2). Populations 1–6 were in the low-hydrothermal-synthesis index region (average 24.97), while populations 7–10 were in the high-hydrothermal-synthesis index region (average 31.98) (Tables 1 and 2). The equation of hydrothermal synthesis index was $S = \displaystyle\sum\limits_{t = 1}^{12} {0.18{r_t}/{{1.045}^t}} $ (r, monthly rainfall; t, monthly mean temperature) (Bailey, 1979). For each population, 17 to 30 individuals (spaced 30 m apart) were sampled (Table 3). Young and healthy leaves were randomly sampled from individuals and immediately stored with silica gel in ziplock plastic bags until DNA extractions were carried out. We chose A. lavandulifolia as the out-group.

2.2 Molecular methods

Total genomic DNA was extracted using AxyPrep Genomic DNA Mini Kits (Axygen Inc., Beijing, China) following the manufacturer's instructions. DNA quality was checked on a 1.0% agarose gel. Several pairs of cpDNA primers designed by Hamilton (1999), Taberlet et al. (1991) , and Sang et al. (1997) were used in the initial screening. Two pairs of primers, trnL (5'-CGGAATTGGTAGACGCTACG-3') and trnF (5'-ATTTGAACTGGTGACACGAG-3') (Sang et al., 1997 ), identified sequence variations in the sampled individuals and therefore were used for all remaining individuals. Polymerase chain reaction (PCR) was performed in a 25-μL reaction volume, containing 40 ng of genomic DNA, 1.0 U of Taq polymerase (Axygen Inc., Beijing, China), 3 mmol/L MgCl₂, 500 μmol/L each dNTP, 20 mmol/L Tris-HCl (pH 8.3), 100 mmol/L KCl, and 0.3 μmol/L primer. The amplification condition was an initial denaturation step at 94 °C for 3 min, followed by 30 cycles of 30 s at 94 °C, 30 s at 55 °C, 1 min at 72 °C, and a final 5-min extension step at 72 °C. The PCR products were determined by 1.0% agarose gel electrophoresis. The amplification products were purified using an AxyPrep PCR Purification Kit, following the manufacturer's protocol (Axygen Inc., Beijing, China). Purified DNA was sequenced by the MEIJI sequencing company in Shanghai, China, applying the PCR-primers as sequencing primers.

2.3 Data analyses

DNA sequences were aligned using the CLUSTAL X program (Thompson et al., 1997 ), with subsequent manual adjustments in MEGA4 (Tamura et al., 2007 ). A matrix of combined sequences was constructed for the 243 individuals that we examined, and different cpDNA sequences were identified as haplotypes.

Basic population genetic parameters were estimated for three groups of populations: the low-hydrothermal-level region (populations 1–6); the high-hydrothermal-level region (populations 7–10); and finally, all populations. All parameters were calculated with D_NASP 5.10.01 (Librado and Rozas, 2009), including the number of segregating sites (S), the number of haplotypes (N_h), the haplotype diversity (H_d), the average number of nucleotide differences per site between two sequences in a sample, π (Nei and Li, 1979; Nei, 1987), and the average number of pairwise nucleotide differences (k).

Phylogenetic analyses of cpDNA haplotypes were performed with maximum parsimony (MP), using PAUP version 4.0 (Swofford, 2002). Heuristic search was implemented with 100 random additional sequence replicates, tree-bisection-reconnection (TBR) branch swapping, MULPARS option, and ACCTRAN optimization. To evaluate the relative robustness of the clades found in the most parsimonious tree, bootstrap analysis was conducted using 1,000 replicates with a simple taxon addition. Genetic differentiation among populations at the three different sampling levels was estimated by pairwise F_ST values (Wright, 1951). AMOVA was performed to analyze the source of variation among populations, using Arlequin 3.0 (Excoffier et al., 2005 ) with 1,000 replicates of bootstrap.

3 Results 3.1 Haplotype distribution and genetic diversity

Sequence data were obtained for one loci from on average 101, 142, and 243 individuals from the low-hydrothermal-level region, the high-hydrothermal-level region, and the entire species' range, respectively. The length of the aligned trnL–trnF DNA sequences (including trnL and the trnF spacer region) ranged between 849 and 863 bp with two insertions. The analysis of cpDNA variation identified seven haplotypes (HapA-HapG) (Table 4). Haplotype C was the most abundant, occurring in three populations, followed by haplotypes A, B, D, and G, which occurred in two populations; and the remaining haplotypes were found in only a single population (Table 3 and Figure 1). We confirmed the division of the native range into two areas using hydrothermal synthesis index data (Table 3). Six populations in the low-hydrothermal-level region of the species were dominated by four different haplotypes. The other four populations in the high-hydrothermal-level region of the species were dominated by three different haplotypes (Table 3 and Figure 1).

We identified a total of two, three, and three segregating sites for the low-hydrothermal-level region, the high-hydrothermal-level region, and all populations, respectively. Summary statistics of sequence variation are given in Table 5. Overall haplotype diversity (H_d) and nucleotide sequence (π) diversity for A. halodendron were 0.706±0.001 and 0.0013±0.0001, respectively. At the regional level, haplotype diversity and nucleotide diversity between two regions varied between 0.318 (low-hydrothermal-level region) and 0.671 (high-hydrothermal-level region), and between 0.0006 (low-hydrothermal-level region) and 0.0015 (high-hydrothermal-level region), respectively. The population from the high-hydrothermal-level region had higher haplotype diversity and nucleotide diversity than the population from the low-hydrothermal-level region.

3.2 Phylogenetic analyses and genetic structure

Maximum parsimony analysis resulted in a single tree (length=14 and consistency index =0.5834). Two clades with low bootstrap support were identified: one consisting of the Hap C, G, F, E, and D, with four haplotypes (Hap C, G, F, and Hap D) distributed in the low-hydrothermal-level region and one haplotype (Hap E) in the high-hydrothermal-level region; and the other consisting of the rest of Hap A and Hap B, both occurring in the high-hydrothermal-level region (Figure 2).

The analysis of molecular variance (AMOVA) showed that around 17.54% of the variation was attributed to the low-/high-hydrothermal-level regions differentiation, and between- populations variation accounted for just over one-half (55%) of the total variation, indicating very strong differentiation and little gene flow between regions and populations (Table 6). Within the low- and high-hydrothermal-level regions, between-populations differentiation was respectively 81.52% and 93.77% (Table 6). Overall, these results strongly indicate that haplotypes are geographically structured across the species' distribution range. F_ST values were 0.728, 0.815, and 0.938 at the species level, low-hydrothermal-level region, and high-hydrothermal-level region, respectively.

4 Discussion and conclusion

In this study, we examined nucleotide variation at two nuclear loci of A. halodendron sampled from the Horqin Sandy Land. We confirmed the division of the native range into two areas, using hydrothermal synthesis index data (Table 3). The analysis of cpDNA variation identified seven haplotypes (Table 4 and Figure 2). Six populations in the low-hydrothermal-level region of the species were dominated by four different haplotypes. The other four populations in the high-hydrothermal-level region of the species were dominated by three different haplotypes (Table 3 and Figure 2).

Many studies have demonstrated that endemic species tend to possess high levels of genetic diversity (Wang et al., 2010 ; Ge et al., 2011 ). Compared with similar studies (Wang et al., 2011 ), the level of haplotype diversity (H_d=0.706) and the nucleotide diversity (π=0.0013) within the whole populations of A. halodendron were somewhat lower than that of the endemic plants in northern of China. Such a low H_d value and π value indicated that this species was very adapted to the sandy land environment, and this character might have made it become the dominant and constructive species in Horqin Sandy Land. At the same time, the populations from the low-hydrothermal-level region exhibited lower haplotype diversity and more limited nucleotide diversity, and were a subset of that observed in the populations from the the high-hydrothermal-level region. These results lend support to the scenario described by Huang et al. (2011) using ISSR markers and Huang et al. (2013) using psbA-trnH. This fact might be due to environmental differences (Li et al., 2009 ; Bao et al., 2010 ; Liu et al., 2010 ; Li and Wu, 2011). The low-hydrothermal-level region is a region of severe desertification (Zhao et al., 2000 ; Bao et al., 2010 ). In the process of ecological restoration, species increasing adapt to the local environment, which, in turn, reduces the genetic diversity (Wang et al., 2009 ; Huang et al., 2011 ).

Phylogenetic relationships analysis of the cpDNA sequences collected from A. halodendron in the different hydrothermal-level regions in Horqin Sandy Land did not entirely cluster according to populations or hydrothermal-level regions. This result is similar to that of some other studies (Chen et al., 2009 ; Hu et al., 2010 ; Zhou et al., 2010 ; Zhang et al., 2015 ). The MP tree indicated that the seven haplotypes formed two clades: the haplotypes in clade I all came from the high-hydrothermal-level region (63% bootstrap support); the haplotypes (except Hap E) in clade II came from the low-hydrothermal-level region (58% bootstrap support). This cluster process was closely related to the hydrothermal gradients in Horqin Sandy Land. The results of phylogenetic relationships analysis show that hydrothermal conditions change determines the increase of the haplotypes distribution difference.

The spatial analyses of genetic variation in A. halodendron in the different hydrothermal-level regions in Horqin Sandy Land indicate that the between-populations genetic differentiations are high within the low-, the high-hydrothermal-level region, and the total distribution range (Table 6). Our data revealed significant genetic differences among populations in different hydrothermal level regions. The increase in population genetic differentiation might be related to strong human disturbance, high rates of habitat fragmentation, and decreasing population size in Horqin Sandy Land (Zhao et al., 2000 ; Wang et al., 2010 ). It is generally accepted that genetic differentiation would have increased among the native populations due to the changed environment (Wang et al., 2009 ).

In conclusion, there are significant genetic differences between A. halodendron populations from different hydrothermal-level regions. This information about genetic variation has important implications for restoring and managing the degraded ecosystems in arid and semi-arid areas. It is particularly important to understand the genetic variation of extant populations over the entire range of the species' distribution. Additional research is needed to determine levels of genetic variation in A. halodendron throughout its entire distributional range.

Acknowledgments:

The authors thank all the members of Naiman Desertification Research Station, China Academy of Sciences (CAS), for their help in the field work. We acknowledge the China Meteorological Administration (Beijing, China) for help on the meteorological data information support. This study was financially supported by research projects 2016YFC0500907, 2017FY100205, 41201561, Y551821001, and 145RJYA269.