Population Demographic History of a Rare and Endangered Tree Magnolia sprengeri Pamp. in East Asia Revealed by Molecular Data and Ecological Niche Analysis

Quaternary climate and environment oscillations have profoundly shaped the population dynamic history and geographic distributions of current plants. However, how the endangered and rare tree species respond to the climatic and environmental fluctuations in the subtropical regions of China in East Asia still needs elucidation. In this study, we collected 36 natural populations of an endangered and rare tree species Magnolia sprengeri Pamp. in subtropical China to determine the demographic history, and modeled the changes of geographic distributions of this species in East Asia based on the MaxEnt ecological niche analyses. In addition, we sequenced three maternally inherited chloroplast DNA fragments (matK, trnH-psbA, and rbcL) for all the natural populations which covered the whole geographic distributions of M. sprengeri. Population genetic analysis showed that the endangered tree species have a low level of chloroplast DNA diversity. However, the genetic variation contribution within populations was greater than that among populations (FST = 0.276), which demonstrated a high level of genetic differentiation. Interestingly, some unique chloroplast DNA haplotypes and higher genetic variations were identified in the Qinling-Daba Mountains, Central China, and Tianmu Mountains of Zhejiang province, East of China in East Asia. Combining with the species distribution modeling, we speculated that these areas might be the potential glacial refugia for the endangered plant M. sprengeri. Phylogeographic analysis demonstrated that the geographic factors (e.g., mountains, rivers, and other isolation barriers) had little effect on the genetic divergence among populations. Ecological niche modeling further revealed that the natural populations of M. sprengeri did not experience significant geographic distribution changes from the last glacial maximum to the present time. These findings are in line with the analysis results of the multimodal mismatch patterns of the chloroplast DNA variations. To protect the endangered species M. sprengeri, in situ and ex situ conservation strategies should be formulated for the natural populations with higher genetic variations.


Introduction
During the Quaternary periods, the global climate experienced repeated cycles of glacial and interglacial stages [1], which lead to the changes of geographic distributions Recently, some chloroplast genes were widely used to detect the population history of Magnoliaceae species [31,32]. For instance, Massoni et al. reconstructed the phylogenetic relationships of Magnoliaceae with 12 molecular markers including matK and rbcL from plastid genome [33]. Some other studies used the chloroplast intergenic regions (e.g., psbA-trnH) to determine the evolutionary relationships of the family Magnoliaceae [34,35]. They found that these chloroplast DNA fragments (i.e., matK, trnH-psbA and rbcL) can be easily obtained by PCR amplification and sequencing for Magnoliaceae species [34,35]. In this study, we investigated the demographic history and the species distribution fluctuations of M. sprengeri using cpDNA sequence data (matK, trnH-psbA and rbcL) and ecological niche modeling. The main objectives of this study are: (I) to examine the population dynamic history of M. sprengeri; (II) to determine the species distribution changes of M. sprengeri during the different historical periods based on the MaxEnt ecological analyses; and (III) to provide a scientific basis for the protection and management of the long-lived endangered species.

Field Investigation and Samples Collection
During the field investigations, we sampled 36 natural populations of M. sprengeri, which covered the whole geographic distributions in subtropical regions of China in East Asia. However, due to the climatic fluctuations and recent human activities, the samples for each population were less than merely 87 M. sprengeri individuals were collected from the eight provinces (Chongqing, Gansu, Guizhou, Hubei, Hunan, Shaanxi, Sichuan, Zhejiang) in China. Healthy and tender leaves were collected in the field and immediately put into a sealed bag containing silica gel for drying. At least 100 m interval between individuals was taken during sampling, and the altitude, longitude, and latitude information of each sampling site was recorded in detail. The vouchers of all materials and documents were deposited in the College of Life Sciences, Northwest University, Xi'an, China (Table S1).

Genetic Diversity and Population Structure Analyses
All DNA sequences were aligned using CLUSTAL X [41] and edited manually in BioEdit v7.2.5 [42]. We concatenated the three cpDNA alignments for further genetic analyses according to their maternally inherited characteristics. DnaSP v6.12.03 [43] was used to calculate basically population genetic parameters, including the number of segregating sites (S), haplotype diversity (Hd), the number of haplotypes (N), and nucleotide diversity (π). Neutrality test statistics of Tajima's D were performed on all samples using DnaSP to assess whether the combined sequences evolved neutrally. We also used mismatch distribution analysis to determine the demographic dynamic changes. Generally, the structures of populations that have experienced demographic expansion are unimodal, while those kept at a stable size exhibit multimodal mismatch distribution structures. In addition, a median-joining network was constructed in NETWORK 10.2.0.0 [44] to infer the relationships of haplotypes. The geographical distributions of cpDNA haplotypes for M. sprengeri were visualized using ArcGIS v10.2 (ESRI, Redlands, CA, USA) in accordance with the methods of Li et al. [45].
A molecular variance analysis (AMOVA) with 1000 permutations was conducted using Arlequin v 3.11 [46]. The Mantel test with 999 permutations was conducted using GenALEx 6.502 [47] to identify whether there is a significant association between genetic distance [48] and geographic distance. Using the program PermutCpSSR_1.2.1 [49], the population differentiation was measured by G ST (gene differentiation coefficient) and N ST (genetic variation coefficient influenced by both the haplotype frequencies and genetic distances between haplotypes), with respect to the haplotypes, and compared by a test with 1000 permutations.
Phylogenetic analyses of the identified haplotypes were performed to detect the evolutionary relationships of chloroplast haplotypes. Maximum-likelihood (ML) tree generation and bootstrap analyses were performed using the program RAxML [50]. We selected the best-scoring ML tree using a generalized time reversible plus gamma model of sequence evolution with 1000 bootstrap replicates. Sequences of the same cpDNA region of Magnolia liliiflora Desr. were used as the outgroup according to previous phylogenetic results [51]. To obtain the Bayesian inference cladogram, we constructed the phylogenetic tree of haplotypes in MRBayes 3.2.2 [52] based on GTR+I model and retained every 300 generations from 3,000,000 random tree rotations. The results were visualized using Figtree v.1.3.1 [53].

Species Distribution Modeling
Data on bioclimatic environmental variables were downloaded from the WorldClim Website (http://www.worldclim.org/ accessed on 5 January 2021). The climatic data for each period included 19 bioclimatic variables (bio1-19) with a resolution of 2.5 arc-minute. The future climate data used in this study (the year 2050 and 2070) are based on the CCSM4 model with strong simulation capability in China [54], which includes four emission scenarios of RCP (Representative Concentration Pathways) 2.6, RCP4.5, RCP6.0 and RCP8.5 in the fifth IPCC (Intergovernmental Panel on Climate Change) emission report [55].
By searching China National Specimen Resource Platform (NSII, http://www.nsii. org.cn// accessed on 5 January 2021) and relevant literature to determine the geographic data of M. sprengeri, combing with the field investigation records, a total of 98 record points of M. sprengeri populations were obtained after screening the samples for authenticity and excluding over-dense loci.
The whole 19 bioclimatic variables were used to simulate the climate factors with a large contribution. These operations relied on MaxEnt_version 3.4.1 software [56]. Modeling of the present distribution of M. sprengeri was undertaken using the maximum entropy algorithm, with default parameter settings (maximum number of background points 10,000) [57]. Measurable environmental factors which affect M. sprengeri's distribution were calculated by creating response curve and doing jackknife, following the methodology of Cao et al. [58]. Statistical Package for Social Sciences, version 24 (IBM Corp, Armonk, NY, USA) was used for correlation analysis of environmental variable information. Subsequently, environment variables whose correlation is greater than 0.8 and less than −0.8 were removed in order to avoid the over-fitting phenomenon caused by the high correlation of environment variables. Combined with the variables whose correlation values are much larger than the mean value, the following eight variables (bio2, 3,4,9,10,11,12,15) are finally obtained.
In a previous study, all of the ecological niche models simulated by MaxEnt had high predictive capacities. In addition, these models were generally good predictors of species' The area under the ROC curve (AUC) for a random classifier is 0.5, while that for a perfect classifier is 1 [59]. Through Conversion Tools-ASCII to Raster in ArcGIS, the simulation results obtained after MaxEnt running optimized parameters were converted into Raster data and the suitable area was divided. MaxEnt was used to predict the distribution of M. sprengeri in the contemporary, Last Glacial Maximum (LGM) and Last Interglacial Age (LIG).

The Distribution and Relationships of Chloroplast DNA (cpDNA) Haplotypes
In this study, three cpDNAs fragments (matK, trnH-psbA and rbcL) from 87 M. sprengeri individuals (belonging to 36 populations) were aligned ( Figure 1). The results of three DNA fragments' primers electrophoretic maps showed that the primers could amplify the target bands in M. sprengeri ( Figures S1-S3). The total alignment length was 1650 bp, and 10 nucleotide substitutions revealed six haplotypes (HP1-HP6) (Table S3, Figure S4). For the whole population, a low level of haplotype diversity (Hd = 0.196) was detected. Almost all populations contain only one haplotype while the populations 7 and 15 had multiple haplotypes. Haplotype HP1 was the most common, found in 34 populations, while endemic haplotypes HP3 and HP4 were unique haplotypes and only distributed in Hubei province, central China. Haplotype HP1 may be an ancestral chloroplast haplotype, giving rise to other haplotypes due to its location in the central position of network analyses. Haplotypes HP2 and HP3 are congregated on a stem in the NETWORK result. Interestingly, three populations 8, 9, and 17 possess the unique cpDNA haplotypes, which are distributed in the boundary between southern Hubei and Hunan in central China. and less than −0.8 were removed in order to avoid the over-fitting phenomenon caused by the high correlation of environment variables. Combined with the variables whose correlation values are much larger than the mean value, the following eight variables (bio2, 3,4,9,10,11,12,15) are finally obtained.
In a previous study, all of the ecological niche models simulated by MaxEnt had high predictive capacities. In addition, these models were generally good predictors of species' occurrences according to the area under the receiver-operating characteristic (ROC) curve. The area under the ROC curve (AUC) for a random classifier is 0.5, while that for a perfect classifier is 1 [59]. Through Conversion Tools-ASCII to Raster in ArcGIS, the simulation results obtained after MaxEnt running optimized parameters were converted into Raster data and the suitable area was divided. MaxEnt was used to predict the distribution of M. sprengeri in the contemporary, Last Glacial Maximum (LGM) and Last Interglacial Age (LIG).

The Distribution and Relationships of Chloroplast DNA (cpDNA) Haplotypes
In this study, three cpDNAs fragments (matK, trnH-psbA and rbcL) from 87 M. sprengeri individuals (belonging to 36 populations) were aligned ( Figure 1). The results of three DNA fragments' primers electrophoretic maps showed that the primers could amplify the target bands in M. sprengeri ( Figures S1-S3). The total alignment length was 1650 bp, and 10 nucleotide substitutions revealed six haplotypes (HP1-HP6) (Table S3) ( Figure S4). For the whole population, a low level of haplotype diversity (Hd = 0.196) was detected. Almost all populations contain only one haplotype while the populations 7 and 15 had multiple haplotypes. Haplotype HP1 was the most common, found in 34 populations, while endemic haplotypes HP3 and HP4 were unique haplotypes and only distributed in Hubei province, central China. Haplotype HP1 may be an ancestral chloroplast haplotype, giving rise to other haplotypes due to its location in the central position of network analyses. Haplotypes HP2 and HP3 are congregated on a stem in the NETWORK result. Interestingly, three populations 8, 9, and 17 possess the unique cpDNA haplotypes, which are distributed in the boundary between southern Hubei and Hunan in central China.   Figure 2). The tree produced a haplotype phylogenetic relationship similar to the ones produced by the network analysis. Haplotypes HP1, HP4, HP5, and HP6 were clustered into a large genetic lineage with high support value. They had a rather distant relationship with haplotypes HP2 and HP3.
(cpDNA) haplotypes detected in M. sprengeri. The size of the circle in the network graph corresponds to the proportion of haplotype occurrence frequency in all populations. The largest circle represents the haplotype with the most individuals, and the color of the circle matches the haplotype color in the distribution.
The phylogenetic tree was constructed based on ML (maximum likelihood) and Bayesian analysis for the three concatenated cpDNA fragments of M. sprengeri, setting Magnolia liliiflora Desr. (LI) as the outgroup, and the bootstrap values were over 50 at all nodes ( Figure 2). The tree produced a haplotype phylogenetic relationship similar to the ones produced by the network analysis. Haplotypes HP1, HP4, HP5, and HP6 were clustered into a large genetic lineage with high support value. They had a rather distant relationship with haplotypes HP2 and HP3.

Genetic Differentiation and Population Structure
According to genetic variation analysis, we found low levels of haplotype diversity (Hd = 0.196) with the average within-population diversity (HS = 0.068) and the total genetic diversity (HT = 0.169). Additionally, a significantly higher GST than NST (NST = 0.268, GST = 0.595, p < 0.05) was observed, revealing that M. sprengeri did not have a pedigree geographic structure (Table 1).

Genetic Differentiation and Population Structure
According to genetic variation analysis, we found low levels of haplotype diversity (Hd = 0.196) with the average within-population diversity (H S = 0.068) and the total genetic diversity (H T = 0.169). Additionally, a significantly higher G ST than N ST (N ST = 0.268, G ST = 0.595, p < 0.05) was observed, revealing that M. sprengeri did not have a pedigree geographic structure (Table 1). Based on cpDNA variation, mismatch analyses showed the multimodel mismatch distribution, which revealed that M. sprengeri maintained relatively stable population sizes throughout the last glacial period (Figure 3). Arlequin software was used to analyze the molecular variation (AMOVA) of cpDNA sequences ( Table 2). The results showed that genetic variation mainly existed within populations, accounting for 72.39%. F-statistics showed that there was obvious genetic differentiation among populations. Correlation between genetic and geographic distances was detected by the Mantel test (r = 0.247, p = 0.03, 999 permutations), indicating that geographical factors had little effect on the genetic differentiation of M. sprengeri.

Model Performance and Contributions of Variables
All of the four ecological niche models simulated by MaxEnt had high predictive capacities (AUC > 0.95) (Figure 4), and the projected present distribution is consistent with collection records. Some of the changing trends of suitable distributions of M. sprengeri obtained by MIROC (Model for Interdisciplinary Research on Climate) are similar to those obtained by CCSM (Community Climate System Model). Arlequin software was used to analyze the molecular variation (AMOVA) of cpDNA sequences ( Table 2). The results showed that genetic variation mainly existed within populations, accounting for 72.39%. F-statistics showed that there was obvious genetic differentiation among populations. Correlation between genetic and geographic distances was detected by the Mantel test (r = 0.247, p = 0.03, 999 permutations), indicating that geographical factors had little effect on the genetic differentiation of M. sprengeri.

Model Performance and Contributions of Variables
All of the four ecological niche models simulated by MaxEnt had high predictive capacities (AUC > 0.95) (Figure 4), and the projected present distribution is consistent with collection records. Some of the changing trends of suitable distributions of M. sprengeri obtained by MIROC (Model for Interdisciplinary Research on Climate) are similar to those obtained by CCSM (Community Climate System Model).
These eight variables, bio2 (mean diurnal temperature range), 3 (isothermality), 4 (temperature seasonality), 9 (mean temperature of the driest quarter), 10 (mean temperature of the warmest quarter), 11 (mean temperature of the coldest quarter), 12 (annual precipitation), and 15 (precipitation seasonality) showed a higher gain compared to others ( Table 3). The contribution rate and importance of environmental variables to the distribution of M. sprengeri in different periods are shown in Table 4. For every simulated period, the top two environmental factors are bio2 (mean diurnal temperature range) and bio9 (mean temperature of the driest quarter).  These eight variables, bio2 (mean diurnal temperature range), 3 (isothermality), 4 (temperature seasonality), 9 (mean temperature of the driest quarter), 10 (mean temperature of the warmest quarter), 11 (mean temperature of the coldest quarter), 12 (annual precipitation), and 15 (precipitation seasonality) showed a higher gain compared to others ( Table 3). The contribution rate and importance of environmental variables to the distribution of M. sprengeri in different periods are shown in Table 4 (Table 4). For every simulated period, the top two environmental factors are bio2 (mean diurnal temperature range) and bio9 (mean temperature of the driest quarter).    The results of the jackknife test of variables' contribution are shown in Figure 5. Bio9 (mean temperature of the driest quarter) and bio11 (mean temperature of the coldest quarter) provided very high gains (>1.8) when used independently, indicating that Bio9 and bio11 contained more useful information by themselves than the other variables did. However, bio2, bio4, bio12, and bio15 had moderate gain when used independently. Other variables, including bio3 and bio10, had low gains when used in isolation; they did not contain much information by themselves. The results of the jackknife test of variables' contribution are shown in Figure 5 ( Figure 5). Bio9 (mean temperature of the driest quarter) and bio11 (mean temperature of the coldest quarter) provided very high gains (>1.8) when used independently, indicating that Bio9 and bio11 contained more useful information by themselves than the other variables did. However, bio2, bio4, bio12, and bio15 had moderate gain when used independently. Other variables, including bio3 and bio10, had low gains when used in isolation; they did not contain much information by themselves. The relationship between the existence probability of M. sprengeri and environmental factors was investigated according to the response curve of environmental factor variables ( Figure 6). When the existence probability of M. sprengeri is greater than 0.5, the corresponding environmental factor value is beneficial to the growth of M. sprengeri. The relationship between the existence probability of M. sprengeri and environmental factors was investigated according to the response curve of environmental factor variables ( Figure 6). When the existence probability of M. sprengeri is greater than 0.5, the corresponding environmental factor value is beneficial to the growth of M. sprengeri.

M. sprengeri Potential Distribution
Using the natural breaks method, the potential distribution of M. sprengeri was divided into four grades (not suitable, marginally suitable, moderately suitable, and highly suitable areas) (Table 5, Figure 7). Our simulation showed that the potential species distribution range was continuous from 20 • N to 35 • N during the LIG (last interglacial). During the LGM (Last Glacial Maximum), when temperature decreased, our results showed that M. sprengeri had the widest distribution and all of the suitable areas are within the subtropical areas of China, with highly suitable areas in the Qinling-Daba Mountains and Sichuan Basin. The projection of the model over the present bioclimatic conditions showed that the habitat was consistently suitable in subtropical China (ca. 22-34 • N) for M. sprengeri, although some occurrences were outside the predicted distribution with marginally and moderately suitable probability. At present, the area of highly suitable and moderately suitable areas is reduced, and the highly suitable area at the Sichuan Basin has almost disappeared, while the proportion of the marginally suitable area increases. From LIG to the present, the distribution area of M. sprengeri has moved northward and expanded during the glacial climate. 021, 12, x FOR PEER REVIEW 10 of 18

M. sprengeri Potential Distribution
Using the natural breaks method, the potential distribution of M. sprengeri was divided into four grades (not suitable, marginally suitable, moderately suitable, and highly suitable areas) ( Table 5) (Figure 7). Our simulation showed that the potential species distribution range was continuous from 20° N to 35° N during the LIG (last interglacial). During the LGM (Last Glacial Maximum), when temperature decreased, our results showed that M. sprengeri had the widest distribution and all of the suitable areas are within the subtropical areas of China, with highly suitable areas in the Qinling-Daba Mountains and Sichuan Basin. The projection of the model over the present bioclimatic conditions showed that the habitat was consistently suitable in subtropical China (ca. 22-34° N) for M. sprengeri, although some occurrences were outside the predicted distribution with marginally and moderately suitable probability. At present, the area of highly suitable and moderately suitable areas is reduced, and the highly suitable area at the Sichuan Basin has almost disappeared, while the proportion of the marginally suitable area increases. From LIG to the present, the distribution area of M. sprengeri has moved northward and expanded during the glacial climate.

Suitable Distribution under Future Climate Scenarios
The MaxEnt model is used to predict the potential geographical distribution changes of M. sprengeri under four emission scenarios of RCP2.6, RCP4.5, RCP6.0, and RCP8.5 in the future (2050 and 2070) ( Table 6). The potential distribution of M. sprengeri shows a shrinking trend in the future climatic environment (Figure 8). Compared with the current potential distribution area, the distribution area of this species in eastern China will be greatly reduced by 2050, gradually approaching the Qinling-Daba Mountains and Dabie Mountains. With the prediction of future climate change scenarios changing from optimistic to pessimistic, the potential future distribution area of M. sprengeri shrinks to a greater extent. At the same time, from 2050 to 2070, the potential future distribution area of M. sprengeri decreases greatly. The potential distribution of M. sprengeri was divided into four grades by the natural breaks method. Gray, green, yellow, and red areas represent not suitable, marginally suitable, moderately suitable, and highly suitable areas, respectively.

Suitable Distribution under Future Climate Scenarios
The MaxEnt model is used to predict the potential geographical distribution changes of M. sprengeri under four emission scenarios of RCP2.6, RCP4.5, RCP6.0, and RCP8.5 in the future (2050 and 2070) ( Table 6). The potential distribution of M. sprengeri shows a shrinking trend in the future climatic environment (Figure 8). Compared with the current potential distribution area, the distribution area of this species in eastern China will be greatly reduced by 2050, gradually approaching the Qinling-Daba Mountains and Dabie Mountains. With the prediction of future climate change scenarios changing from optimistic to pessimistic, the potential future distribution area of M. sprengeri shrinks to a greater extent. At the same time, from 2050 to 2070, the potential future distribution area of M. sprengeri decreases greatly. Table 6. Four emission scenarios using in this study.

Emission
Description The radiative forcing rose to 8.5 W/m 2 , and the CO 2 equivalent concentration reached about 1370 mL/m 3 in 2100.

RCP6.0
The radiative forcing stabilized at 6.0 W/m 2 , and the CO 2 equivalent concentration stabilized at about 850 mL/m 3 after 2100.

RCP4.5
The radiative forcing stabilized at 4.5 W/m 2 , and the CO 2 equivalent concentration stabilized at about 600 mL/m 3 after 2100.

RCP2.6
The radiative forcing reached its peak before 2100 and decreased to 2.6 W/m 2 by 2100. The peak CO 2 equivalent concentration was about 490 mL/m 3 .

RCP6.0
The radiative forcing stabilized at 6.0 W/m , and the CO2 equivalent concentration stabilized at about 850 mL/m 3 after 2100.

RCP4.5
The radiative forcing stabilized at 4.5 W/m 2 , and the CO2 equivalent concentration stabilized at about 600 mL/m 3 after 2100.

RCP2.6
The radiative forcing reached its peak before 2100 and decreased to 2.6 W/m 2 by 2100. The peak CO2 equivalent concentration was about 490 mL/m 3 .
1 Representative Concentration Pathways.  ). The green area represents the overlap of the present and the predicted ranges; blue represents potential range expansion; and red represents potential range contraction.

Diversity and Genetic Structure
Genetic diversity is an important component of biodiversity, which has important ecological influence on natural populations [60]. The higher the genetic diversity or the richer the genetic variation of a species, the stronger ability it has to adapt to environmental changes, and the easier it is to expand its distribution ranges and develop new environments [61]. In this study, we found a generally low level of genetic diversity of the endangered tree M. sprengeri based on cpDNA datasets, which is similar to the study results of other endangered plants Dunnia sinensis in East Asia [62]. Commonly, the genetic diversity of a species is influenced by many factors, including evolutionary history, geographical distribution, and biological characteristics of the species itself [63]. Among the nine geographic regions of M. sprengeri in this study, the genetic diversity of M. sprengeri in Hubei and Hunan populations was the highest, followed by Guizhou province. Therefore, it was inferred that M. sprengeri gradually spread from the border area of Hubei Province to the south. Since the establishment of a new population, the gene frequency of the new population has depended on the genotype of the first few or dozens of individuals. When encountering the influence of adverse external conditions along the way, the gene frequency of the offspring will change with the change of gene frequency of the surviving individual. The current results showed that the boundary region of Hubei and Hunan Provinces and the earlier invaded region of Guizhou Province had higher genetic diversity, while the newly invaded region had lower genetic variations. Additionally, by selecting samples for distribution simulation, M. sprengeri is found to be mainly distributed in the subtropical monsoon climate area [64]. The ecological analysis  (d,h). The green area represents the overlap of the present and the predicted ranges; blue represents potential range expansion; and red represents potential range contraction.

Diversity and Genetic Structure
Genetic diversity is an important component of biodiversity, which has important ecological influence on natural populations [60]. The higher the genetic diversity or the richer the genetic variation of a species, the stronger ability it has to adapt to environmental changes, and the easier it is to expand its distribution ranges and develop new environments [61]. In this study, we found a generally low level of genetic diversity of the endangered tree M. sprengeri based on cpDNA datasets, which is similar to the study results of other endangered plants Dunnia sinensis in East Asia [62]. Commonly, the genetic diversity of a species is influenced by many factors, including evolutionary history, geographical distribution, and biological characteristics of the species itself [63]. Among the nine geographic regions of M. sprengeri in this study, the genetic diversity of M. sprengeri in Hubei and Hunan populations was the highest, followed by Guizhou province. Therefore, it was inferred that M. sprengeri gradually spread from the border area of Hubei Province to the south. Since the establishment of a new population, the gene frequency of the new population has depended on the genotype of the first few or dozens of individuals. When encountering the influence of adverse external conditions along the way, the gene frequency of the offspring will change with the change of gene frequency of the surviving individual. The current results showed that the boundary region of Hubei and Hunan Provinces and the earlier invaded region of Guizhou Province had higher genetic diversity, while the newly invaded region had lower genetic variations. Additionally, by selecting samples for distribution simulation, M. sprengeri is found to be mainly distributed in the subtropical monsoon climate area [64]. The ecological analysis showed that temperature difference was the main geographical and climatic factor affecting the genetic variation of M. sprengeri.
Genetic structure represents the distribution patterns of genetic variation within and among populations, which is mainly affected by internal and external factors. The internal factors mainly include gene flow, genetic drift, bottleneck effect, breeding methods, etc., while the external factors mainly include geological history changes, human activity and excessive mining, etc. [65,66]. The results of AMOVA analysis showed that the genetic variation of M. sprengeri mainly presented within the populations, accounting for 72.39% of total variations, and there was a significant genetic differentiation between populations (F ST = 0.276), indicating less gene communication among different geographical populations but more frequent within the populations [67]. The Mantel tests show that geographical factors have little influence on the genetic differentiation of M. sprengeri. Some previous studies found that the seed dispersal mechanism has a significant impact on genetic differentiation among natural populations [68]. The seeds of M. sprengeri are thicker; therefore, the water permeability is poor, often leading to the low seed germination rate, which was also observed in the study of Magnolia dealbata [69]. Plus, seeds of Magnoliaceae plants are red; thus, the seeds that fall to the ground are conspicuous and frequently eaten by pests or animals [70]. Its upgrading ability was also seriously insufficient due to the lack of updated seedlings of the internal basic population. These mechanisms might have caused significant genetic differentiation [71]. Meanwhile, He et al. speculated that the birds and animals living on the fruits of wild cherries may be an important factor affecting the genetic structure of cherries [72]. Additionally, M. sprengeri, as a substitute for precious Chinese herbal medicine [73], has been cut down in large numbers by the long-term man-made destruction, which results in many populations in single plant distribution or fragmented and has further caused the high genetic differentiation and low genetic variations [74].

Population History and Species Distribution Fluctuations
The geographical distribution of species populations is closely related to the environmental impact of the Quaternary Glacial Age [75]. During the ice age, species may have located in refuges in areas that did not experience ice ages, and thus had longer evolutionary history and greater genetic diversity. In the present study, MaxEnt was used to predict the distribution of M. sprengeri in the contemporary, LIG, and LGM in this study. Since many parameters in MIROC are more suitable than CCSM for the simulation of East Asia [76], we referred more to the former model. The results showed that M. sprengeri was the most widely distributed in LGM, which indicated that the population expansion event occurred in the cold climate period of the Last Glacial Maximum. The optimal distribution area existed in the Qinling-Daba Mountains (QDM), the southeast of the margin of Sichuan Basin, and part of the west of the margin of Sichuan Basin. In addition, this species is also distributed in the Sichuan Basin, east of the North China Plain and the middle and lower reaches of the Yangtze River basin. Many northerly distribution tree species are known to be habitat generalists which are able to survive in various habitat types in temperate regions [77]. Combining the forecast area of each period, it was found that the suitable growth altitude of M. sprengeri is 1300-2400 m [25], which is greatly affected by temperature and climate factors. We speculated that M. sprengeri with considerable tolerance has survived the ice age with a slightly total distribution decrease.
Additionally, the occurrence of areas or centers of endemism is commonly attributed to the existence of suitable refugia [65]. The QDM is located in central China separated by the southern subtropical and northern temperate regions, with complex topography, climate, and ecological diversity [66]. In this study, populations (number 7 and 15) with multiple haplotypes and populations (number 8, 9 and 17) with endemic haplotypes (HP3, 4) were distributed in the southeastern edge of the Sichuan Basin (the junction area of Hunan-Hubei provinces and Guizhou Province). It is speculated that the QDM and the whole southeastern edge of the Sichuan Basin were the glacial refugees of M. sprengeri during LGM. However, the predicted highly suitable areas of this species are now decreasing in fragmentation (138 to 116 km 2 , from LGM to the present), possibly due to the environmental changes and recently human activities [78,79]. The predicted highly suitable areas in the Sichuan Basin have been disappeared currently, and the total suitable areas were significantly lower than that of in LGM. Recent studies showed that the distribution area of Euscaphis japonica in East Asia has gradually shrunk as a result of the regional influence of climate change, together with deforestation [80,81]. In this study, a similar result was found that the decrease of the current distribution area of M. sprengeri might be caused by the frequent human activities, the instability of climate and ecology, and the habitat destruction by human beings for economic benefits. Researchers are now realizing that climate change may be a major threat to biodiversity in the next 100 years [82]. Therefore, M. sprengeri's survival situation will gradually improve with the strengthening of human environmental protection consciousness and the earth's ecological restoration.

Conservation Strategies and Implications
The floras with high levels of endemism are relatively vulnerable, more likely to face the risk of habitat losses and extinction [83], indicating that conservation measures for endangered plants are a matter of urgency [84]. Liu et al. conducted a systematic investigation of Magnoliaceae plants in 14 provinces and regions of China and found that the distribution area of most Magnoliaceae plants shrunk with the destruction of habitat, and M. sprengeri is no exception [85]. It also has been listed as a vulnerable endangered species in the Red List [30]. For a long time, M. sprengeri has been used as high-quality wood and medicine [86], and its natural population has been greatly destroyed [87]. The geography of the pedigree of M. sprengeri in research results in differentiation within populations being much higher than that between populations. Therefore, we recommend local protection of the wild natural population of this species [88], especially in the natural population of potential glacial refugia (populations in QDM and the southeastern edge of Sichuan Basin), in order to preserve their high genetic diversity. In addition, it is advised to focus on the protection of groups with endemic haplotypes (populations number 7,8,9,15,17) to prevent loss of M. sprengeri within the population diversity and the risk of extinction. Moreover, populations in different regions can also be protected reasonably and exoterically to increase the gene exchange between populations, improve their genetic diversity, and enhance the survivability of M. sprengeri against adverse environments [89]. Chen et al. suggested that effective gene flow can be spread through seeds [90]. Therefore, the mature seeds of each population can be collected and artificially sown to other populations by employing the cross-sowing method [91], so as to improve the fragmentation of habitat, strengthen gene communication among populations, and improve the level of genetic diversity of wild populations.
As one of the most primitive taxa of primitive angiosperms, the ancestors of M. sprengeri had an abundant genetic basis and genetic variation [92]. However, the rare and endangered species often have less genetic diversity than that of other widespread species of the family Magnoliaceae [93,94], indicating that the process of evolution may not be the main cause of M. sprengeri's genetic diversity. The relatively low genetic diversity may be related to the distribution decrease and habitat fragmentation. In addition, Germplasm resources and artificial seedling technology of M. sprengeri should be established [95] to improve the natural regeneration capacity of M. sprengeri and expand the range of its cultivation. Ex-situ conservation should be carried out for some populations with low genetic diversity [96], and the genetic materials of different populations should be collected as far as possible to promote the extensive exchange of genetic information among different populations.

Conclusions
In this study, population genetic analysis showed that the endangered and rare tree Magnolia sprengeri have a low level of chloroplast DNA diversity and high genetic differentiation among populations. Ecological niche modeling demonstrated that the natural populations of M. sprengeri did not experience significant geographic distribution changes from the Last Glacial Maximum to the present time. The Qinling-Daba Mountains, Central China, and Tianmu Mountains of Zhejiang province, East of China in East Asia might be the potential glacial refugia for the endangered species M. sprengeri. In situ and ex-situ conservation strategies should be formulated for the natural populations of M. sprengeri in the near future.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/f12070931/s1, Figure S1: Electrophoretic map of matK's primers. Figure S2: Electrophoretic map of trnH-psbA's primers. Figure S3: Electrophoretic map of rbcL's primers. Figure S4: All alignments sequences of M. sprengeri used in this study. Table S1: The location and voucher information of M. sprengeri samples used in this study. Table S2: Primers used for DNA barcoding studies. Table S3: Variable sites of the aligned sequences of the three chloroplast DNA (cpDNA) fragments (matK, trnH-psbA and rbcL) from which six cpDNA haplotypes of Magnolia sprengeri were identified.