Population Structure and Genetic Diversity of the Toona ciliata (Meliaceae) Complex Assayed with Chloroplast DNA Markers

Toona ciliata is a deciduous or semi-deciduous tree species and belongs to the Toona genus of the Meliaceae family. Owing to low natural regeneration and over-exploitation, the species is listed as an endangered species at level II in China and its conservation has received increasing concern. Here, we sampled 447 individuals from 29 populations across the range-wide distribution of the T. ciliata complex in China and assessed their genetic variation using two chloroplast DNA markers. The results showed that the overall haplotype diversity and nucleotide diversity per site were high at h = 0.9767 and π = 0.0303 for the psbA-trnH fragment and h= 0.8999 and π = 0.0189 for the trnL-trnL fragment. Phylogenetic analysis supported the division of the natural distribution of T. ciliata complex into western and eastern regions. The genetic diversity was higher in the western region than in the eastern region, showing significant phylogeographic structure. Genetic differentiation among populations was moderate (Φst=42.87%), and the effects of isolation by distance (IBD) were significant. A neutrality test and mismatch distribution analysis indicated that the distribution of the T. ciliata complex generally did not expand, although a few local populations could likely expand after bottleneck effects. The overall results were complementary to and consolidated previous studies using mitochondrial and nuclear DNA markers. We finally discussed strategies for the genetic conservation of the T. ciliata complex.


Introduction
Toona is a genus of deciduous or semi-deciduous broad-leaved trees in the Meliaceae family and is mainly distributed in subtropical and tropical areas [1].The genus includes four species in China, including T. sinensis, T. ciliata, T. macrocarpa, and T. rubriflora, among which T. sinensis and T. ciliata are the most extensively exploited.T. ciliata is widely distributed in southern China, covering Yunnan, Guangxi, Guizhou, Sichuan, Hunan, Hubei, Fujian, and Guangdong Provinces.Based on the phenotypic variation in leaf and floral traits, five varieties of T. ciliata are recognized, including T. ciliata var.ciliata, T. ciliata var.yunnanensis, T. ciliata var.pubescens, T. ciliata var.sublaxiflora, and T. ciliata var.henryi [1].These varieties are geographically distributed in sympatry or parapatry and grow in hills and mountains at elevations between 300 m and 2260 m above sea level.Natural hybridization among varieties cannot be excluded, although a few reports on their hybrids are available in the literature.Because leaf and floral traits are often influenced by environmental factors and exhibit continuous variation within and between species [2,3], the delimitation of varieties of T. ciliata based on these traits remains uncertain.There is still debate about the taxonomic status of these varieties [4].Similar to a previous study that did not specify a variety, here we bring these five varieties together as the T. ciliata complex.Individual samples of T. ciliata were not identified to any specific variety.
T. ciliata is an important fast-growing timber species and is also known as the Chinese mahogany [5][6][7].One remarkable feature is that T. ciliata has a tough texture, straight grain, beautiful pattern, dark reddish-brown heartwood, and lighter sapwood, and hence, it has a significant timber value [8,9].The species also has a potential medicinal value since extracts from its branches and leaves have various biological activities, including antibacterial, antiviral, molluscacide, and antimalarial effects [9].The species is considered as one of the preferred species to adjust the species structure of coniferous and broad-leaved mixed forests in the hills of southern China [9].It is of great significance for developing forest production to vigorously promote the development of T. ciliata.However, due to its low natural regeneration [10] and over-exploitation, T. ciliata is recognized as an endangered species at level II, which is defined as an endangered species of important economic, cultural, or scientific research value.Therefore, its genetic conservation has received increasing concern.
Previous studies on the T. ciliata complex have mainly focused on cultivation, silviculture techniques, ecophysiology, timber properties, and forest community succession [11].For instance, the mixed forests of T. ciliata with Cunninghamia lanceolata, Pinus massoniana, or Pinus yunnanensis could effectively prevent land decline caused by continuous planting, and improve soil fertility [12,13].Genetic studies have focused on provenance trials [14,15], tissue culture [16], and cuttings propagation [17].Population genetic studies with molecular markers include the investigation of genetic diversity and population structure of the T. ciliata complex using SRAP (sequence-related amplified polymorphisms) markers [18], ISSR (inter-simple sequence repeats) markers [19], SSR (simple sequence repeat) markers [20], and nuclear ITS (nuclear internal transcribed spacer) and mitochondrial DNA (mtDNA) markers [4].They showed that there were significant population differentiation and isolation-by-distance (IBD) effects in the natural distribution of the T. ciliata complex in China.Xiao et al. [21] showed that four varieties of T. ciliata (without T. ciliata var.sublaxiflora among the five varieties) in sympatry were genetically well mixed in terms of chloroplast genomes.Other studies using SSR markers include the population genetic structure of T. ciliata var.pubescens [22,23].Analysis of the mating system scored by nuclear SSR markers indicated that the T. ciliata complex is predominantly outcrossing, with partial selfing and inbreeding [24].Recent genome sequencing provides a useful reference for studying the phylogenetic relationship among varieties and for developing nuclear molecular markers for breeding [25,26].All these studies provide a general picture of the genetic diversity and population structure of the T. ciliata complex, which serves as a valuable reference to the genetic conservation of this endangered species in China.
The purpose of this study was to use chloroplast DNA (cpDNA) markers to investigate the population genetic diversity and population structure of the T. ciliata complex in China.This study is complementary to previous studies using mtDNA and unclear SRAP and ITS markers [4,18].It is well known that chloroplast DNA markers have multiple attributes that are suitable for investigating genetic variation within and between species [27,28], including (i) uniparental inheritance (maternal inheritance for angiosperms and paternal inheritance for gymnosperms), (ii) a higher mutation rate per site compared with the plant mitochondrial genome but a lower mutation rate per site compared with the plant nuclear genome, and (iii) a small genome size (120-220 kb) [29].These characteristics ensure that the genetic diversity derived from cpDNA markers is different from that derived from mtDNA and nuclear DNA markers.Moreover, gene flow for the cpDNA markers of angiosperms is mediated by seed flow only, similar to mtDNA markers but different from nuclear markers whose gene flow is mediated by both seed and pollen flow.It is expected in theory that population genetic differentiation derived from cpDNA markers is comparable to that from mtDNA markers if mutation effects are negligible between them [30].Also, it is expected that genetic differentiation derived from cpDNA markers is greater than that from nuclear markers due to smaller genetic drift effects and both seed and pollen flow for the nuclear markers [31].To compare with previous studies [18,21], we investigated range-wide populations of the T. ciliata complex.By synthesizing information on genetic variation from nuclear and organelle genome markers, we assessed the genetic diversity and population structure of the T. ciliata complex from the whole genome perspective.The overall results were then applied in discussing the genetic conservation strategy for the T. ciliata complex.

Sampling and DNA Extraction
The experimental sampling sites covered 11 provinces, and four hundred and fortyseven individuals of 29 populations were collected across the range-wide distribution of the T. ciliata complex in China (Table 1).Here, populations in different localities were determined according to China's administrative division.These sampling sites were the same as those used by Li et al. [18] and Xiao et al. [4] but fewer samples were analyzed.Most leaf samples were collected from 2015 to 2016 [18].Because some leaf samples were used in previous experiments [18], we collected some supplementary samples from provenance trials established at the Zengcheng Experiment Station (113 • 37 ′ E, 23 Approximately 5 g of fresh leaves per sample were immediately dried and preserved on silica gel.Most DNA samples extracted from healthy leaves were provided by Li et al. [18].Supplementary DNA samples were extracted from the leaves of living plants according to the CTAB 2X protocol [32].The quality of DNA extraction was checked by 0.8% (w/v) agarose gel electrophoresis.All quantified DNA samples were stored at −80 • C for subsequent polymerase chain reaction (PCR) amplification.

Primer Screening, Amplification, and Sequencing
According to the literature [33], four pairs of universal cpDNA primers were checked, namely, psbA-trnH, trnT-trnL, trnL-trnL, and trnL-trnF (Table 2).Among them, the intergenic region of psbA-trnH is one of the fastest evolutionary regions in the chloroplast genome, and the sequence has been widely used in genetic diversity analysis [34].After checking the gel bands of the PCR amplifications (Figure S1), we selected two primer pairs, i.e., psbA-trnH and trnL-trnL, for sequencing analyses.PCR amplification was carried out in a 25 µL reaction system, consisting of 1 µL of genomic DNA, 1 µL of forward and reverse primers, 12.5 µL of 2 × ES Taq Mastermix (Yongjin Biotechnology, Guangzhou, China), and 9.5 µL of ddH 2 O.The PCR procedure was carried out in a Dongsheng Thermal Cycler (EDC-810, Suzhou, China).The procedure of the thermal cycling system for both primer pairs was as follows: pre-denaturation for 4 min at 94 • C, followed by 35 cycles of 94 • C for 1 min, 55 • C for 1 min, 72 • C for 1 min, and finally, 72 • C extension for 10 min.The amplification products were stored in the refrigerator at 4 • C and then promptly sent to Sangon Biotech (Shanghai, China) for sequencing.The peak maps obtained from sequencing were viewed and detected by Chromatogram Explorer3.2(Figure S2), which indicated that there were no cluttered peaks or strong base signals.The sequences with high quality and no mixed peak signals were used for downstream analyses.

Analysis of Genetic Diversity
The sequenced fragments from PCR amplifications were aligned using MEGA 7 [35], removing the parts with heterozygous interferences caused by unstable signals at the fronts and ends.Ultimately, we obtained two datasets of sequences for 29 populations, one for the psbA-trnH fragment and the other for the trnL-trnL fragment (see Table S1 in Supplementary Materials for concatenated alignment sequences).We used DnaSP v5 [36] to estimate haplotypes (see Tables S2 and S3 in Supplementary Materials for haplotype sequences), observed haplotype diversity (h), and nucleotide diversity (π) per site for each population and the overall samples.The observed haplotype diversity (h) for each population was estimated by where p i is the frequency of the ith haplotype in the focal population [37].Since many haplotypes were detected, instead of drawing a complex network, we used parsimony tree by MEGA 7 [35] to depict the evolutionary relationships among haplotypes.

Population Genetic Structure
We used DNAsp v5 to estimate genetic differentiation coefficients, including N st [38,39], G st [37], and F st [40].We tested whether N st was larger than G st or not to infer if the phylogeographic structure occurred.Similarly, Arlequin v3.0 [41] was used for molecular analysis of variance (AMOVA) to calculate the distribution of genetic variation within and between populations and to estimate genetic differentiation Φ st [42] between populations.
Isolation by distance (IBD) was tested using regression analysis of (F st /(1 − F st )) on the logarithm of geographic distance [40,43]: A significant difference of the regression coefficient b from zero indicates the presence of IBD effects.Note that the geographical distance between pairwise populations was calculated using their longitude and latitude coordinates (Table 1).A Mantel's test was also conducted to examine the relationship between F st and geographic distance [44].Pearson's correlation (r) between haplotype diversity and elevation was tested for each fragment.
Phylogenetic relationships among 447 individuals of 29 populations were constructed by MEGA 7 [35].The phylogenetic tree among populations was also drawn using the ML (maximum likelihood) method [36].In addition, we investigated the population structure by ADMIXTURE 1.3.0 using the concatenated sequences of the psbA-trnH and trnL-trnL fragments [45].

Population Demography
Tajima's D. [46] and Fu's F. [47] tests were conducted to indicate if a population had experienced expansion after bottleneck effects.Neutrality was tested using Arlequin v3.0 [41] with the two chloroplast DNA fragments.It is expected that both Tajima's D. and Fu's F. are negative if the population has expanded after bottleneck effects.In addition, a mismatch distribution was analyzed for each population.When the observed frequency of pairwise different sites has a single-peaked distribution, consistent with a single-peaked Poisson distribution, population expansion is implied after bottleneck effects.Statistical tests were performed using the sum of squares (SSD) and Harpending's raggedness index (Rag) to check whether the expected SSD (or Rag) was greater than the observed SSD (or Rag).Estimates of other parameters included θ 0 = 2N 1 µ and θ 1 = 2N 2 µ, where N 1 and N 2 are the population sizes before and after population expansion, respectively, and τ(t) is the time elapsed since the sudden expansion of the population.

Haplotype Analysis
The amplified sequences from the psbA-trnH fragment were about 564bp after alignment among 447 individuals.Analysis with DNAsp v5 identified 239 haplotypes among the 447 samples (Table 3).There were 193 haplotypes in total whose abundance was one in overall samples (447 individuals).For simplicity, we used code Hi (i = 1, . .., 46) to represent haplotypes whose abundances were not less than two and code H47 to represent all haplotypes occurring only once in the overall samples.Populations DC, TL, CH, and HD possessed more than 20 haplotypes, while population YF was fixed by haplotype H9.Haplotype H2 was the most frequent haplotype, accounting for 9.4% of the overall samples and occurring in eight populations.Figure 1A shows the evolutionary relationship among the 46 haplotypes.Only a few haplotypes were dominant, and most haplotypes were derived from a few main haplotypes with a couple of mutant bases.
For the sequence of the trnL-trnL fragment, we identified 143 haplotypes from 447 samples (Table 4).There were also 46 haplotypes whose abundances were not less than two and denoted by Hi (i = 1, . .., 46).There were 97 haplotypes that occurred only once in the overall samples.For simplicity, these single-abundant haplotypes were commonly denoted by code H47.Populations XL, WM, and YR possessed more than 20 haplotypes, while populations XJ, SC, NP, JL, and MT were fixed by haplotype H1.Haplotype H1 was the most frequent haplotype, accounting for 29.8% of the overall samples and occurring in 15 populations.Figure 1B shows that most single-abundant haplotypes were derived from a few domain haplotypes with a couple of mutations.#: x + y means that there were x haplotypes whose abundances were not less than two and y haplotypes whose abundance was one in overall samples.
The analysis of haplotype diversity indicated that the observed haplotype diversity ranged from 0 (YF) to 0.9586 (TL) for the psbA-trnH fragment (Table 3).The mean of observed haplotype diversity across populations was 0.7391 ± 0.2477.The nucleotide diversity per site (π) ranged from 0 (YF) to 0.1398 (HS) with a mean of 0.0303 ± 0.0374 (Table 3).The diversity in the overall samples was h = 0.9767 and π = 0.0303.
Figure 2 shows the geographical distribution of haplotypes for the psbA-trnH fragment (A) and the trnL-trnL fragment (B).A generally similar pattern of diversity was observed for both fragments.Populations had larger haplotype diversity in the western region than in the eastern region.An explicit phylogeographic structure was present in the natural distribution of the T. ciliata complex.#: x + y means that there were x haplotypes whose abundances were not less than two and y haplotypes whose abundance was one in overall samples.

Population Genetic Structure
For the psbA-trnH fragment, AMOVA showed that 26.55% of the total genetic variation was from among populations (p-value = 0.0000) and 73.45% was from within populations (Table 5).For the trnL-trnL fragment, AMOVA showed that a major proportion of genetic variation occurred among populations, accounting for 70.02% (p-value = 0.0000).For the concatenated sequences of the psbA-trnH and trnL-trnL fragments, AMOVA showed that 42.84% of the genetic variation occurred among populations (p-value = 0.0000) and 57.16% occurred from within populations.All these results indicated a significant level of difference in the genetic variation among populations of the T. ciliata complex.
A comparison of N st and G st confirmed the presence of phylogeographic structure in the distribution of haplotype diversity.Analysis with the psbA-trnH sequences showed that G st (=0.1546) was significantly smaller than N st (0.2194) (p-value < 5%).Analysis with the trnL-trnL sequences showed that G st (0.2531) was also significantly smaller than N st (0.6972).These results were consistent with the pattern of the haplotype diversity distribution (Figure 2).IBD tests indicated significant but small correlations between genetic differentiation and geographical distances (psbA-trnH: R 2 = 0.0042, p-value = 2.02 × 10 −12 ; trnL-trnL: R 2 = 0.0001, p-value = 9.02 × 10 −10 ) (Figure 3).Pearson's correlation tests indicated that the observed haplotype diversity was not significantly correlated with elevation for the psbA-trnH fragment (r = 0.3528, p-value = 0.0605) but was significantly correlated with elevation for the trnL-trnL fragment (r = 0.4082, p-value = 0.0279).
Figure 4 shows a consensus tree among 29 populations.Generally, the twenty-nine populations could be divided into two groups: the western and eastern regions.The western region covered Yunnan (SM, PW, BS, YR), Guizhou (XY, WM, CH, LD), Guangxi (TL, XL, LL), Sichuan (DC, HD), and Guangdong (YF, LC) Provinces.The eastern region covered Fujian (NP, served for both fragments.Populations had larger haplotype diversity in the western region than in the eastern region.An explicit phylogeographic structure was present in the natural distribution of the T. ciliata complex. (A) (B)

Population Genetic Structure
For the psbA-trnH fragment, AMOVA showed that 26.55% of the total genetic variation was from among populations (p-value = 0.0000) and 73.45% was from within populations (Table 5).For the trnL-trnL fragment, AMOVA showed that a major proportion of genetic variation occurred among populations, accounting for 70.02% (p-value = 0.0000).For the concatenated sequences of the psbA-trnH and trnL-trnL fragments, AMOVA showed that 42.84% of the genetic variation occurred among populations (p-value = Table 5. Analysis of molecular variance (AMOVA) using the psbA-trnH and trnL-trnL sequences and their concatenated sequences of the T. ciliata complex.The analysis of genetic relationships showed a good genetic mixture among the 447 individuals (Figure 5).Many individuals from different populations were more genetically related.A further analysis with ADMIXTURE showed that the optimal number of subpopulations was K = 23, with the minimum cross-validation (CV) error (CV error = 0.1574) (Figure 6).Individuals across populations were genetically well mixed by different proportions of the optimal 23 populations, supporting the phylogenetic relationships among individuals (Figure 5).

Analysis of Population Demography
Neutrality tests with the concatenated sequences of the psbA-trnH and trnL-trnL fragments indicated that most populations had negative Tajima's D or Fu's F values (Table 6).Eight populations had significant Tajima's D values (p-value < 0.05), including populations LL, SM, DC, WM, BS, TL, JX, and PW.Seven populations had significant Fu's F values (p-value < 0.05), including populations DC, YR, LD, TL, CH, HD, and MT.Only two populations (DC and TL) exhibited significant Tajima's D and Fu's F values (negative).There were four populations with positive Tajima's D values, including SC, CB, NP, and JG, suggesting that these populations had not experienced expansion.
Figure 4 shows a consensus tree among 29 populations.Generally, the twenty-nine populations could be divided into two groups: the western and eastern regions.The western region covered Yunnan (SM, PW, BS, YR), Guizhou (XY, WM, CH, LD), Guangxi (TL, XL, LL), Sichuan (DC, HD), and Guangdong (YF, LC) Provinces.The eastern region covered Fujian (NP, WY), Jiangxi (MT, GS, JL, JG), Anhui (HS, JX), Hunan (XN, HP, CB), Zhejiang (SC, XJ), and Hubei (XE) Provinces.Generally, the two regions were connected to each other in the geographical distribution of the T. ciliata complex.The analysis of genetic relationships showed a good genetic mixture among the 447 individuals (Figure 5).Many individuals from different populations were more genetically related.A further analysis with ADMIXTURE showed that the optimal number of subpopulations was K = 23, with the minimum cross-validation (CV) error (CV error = 0.1574) (Figure 6).Individuals across populations were genetically well mixed by different proportions of the optimal 23 populations, supporting the phylogenetic relationships among individuals (Figure 5).

Analysis of Population Demography
Neutrality tests with the concatenated sequences of the psbA-trnH and trnL-trnL fragments indicated that most populations had negative Tajima's D or Fu's F values (Table 6).Eight populations had significant Tajima's D values (p-value < 0.05), including populations LL, SM, DC, WM, BS, TL, JX, and PW.Seven populations had significant Fu's F values (pvalue < 0.05), including populations DC, YR, LD, TL, CH, HD, and MT.Only two populations (DC and TL) exhibited significant Tajima's D and Fu's F values (negative).There were four populations with positive Tajima's D values, including SC, CB, NP, and JG, suggesting that these populations had not experienced expansion.Table 6.The neutrality test and mismatch distribution analysis with the concatenated alignment sequences of the psbA-trnH and trnL-trnL fragments in 29 populations of the T. ciliata complex.Mismatch distribution analyses indicated that populations with significant Tajima's D or Fu's F did not expand significantly, with the probabilities of P (expected SSD ≥ observed SSD) < 0.95 and P (expected RAG ≥ observed RAG) < 0.95 (Table 6).The actual observed curves were essentially multimodal or sawtooth-shaped, indicating that these populations had not recently experienced expansion or a range expansion.Figure 7 shows the mismatch distributions of populations DC and TL with negative Tajima's D and Fu's F values, suggesting that these two populations did not expand after bottleneck effects.Taken together, populations of the T. ciliata complex generally did not experience expansion in its natural distribution in China, although a couple of local populations could exhibit potential expansion.

Genetic Diversity
This study is complementary to previous studies that used both nuclear and mitochondrial DNA markers [4,18,20].The results consolidate previous results of genetic diversity derived from SRAP, SSR, and nuclear ITS markers.Combining the results from organelle (mtDNA, cpDNA) and nuclear genome (SRAP, SSR, and ITS) markers, we conclude two common pa erns regarding the spatial distribution of genetic diversity: (1) The range-wide distribution can be generally divided into western and eastern regions, although the boundaries between them are difficult to delineate.(2) The western region har-

Genetic Diversity
This study is complementary to previous studies that used both nuclear and mitochondrial DNA markers [4,18,20].The results consolidate previous results of genetic diversity derived from SRAP, SSR, and nuclear ITS markers.Combining the results from organelle (mtDNA, cpDNA) and nuclear genome (SRAP, SSR, and ITS) markers, we conclude two common patterns regarding the spatial distribution of genetic diversity: (1) The rangewide distribution can be generally divided into western and eastern regions, although the boundaries between them are difficult to delineate.(2) The western region harbors larger genetic diversity than the eastern region.Note that the western region may be further divided into two regions based on mtDNA markers [4].It is speculated that the T. ciliata complex could likely initialize in the western region and move forward to the eastern side of China, although the original center of T. cilata remains to be determined.This could likely have a similar origination to T. sinensis in China [49].Lu et al. [49] thought that Chinese T. sinensis could originate from the regions of northeastern India and Burma.This needs further verification in future studies.
As expected, cpDNA markers had greater haplotype diversity and nucleotide diversity per site than mtDNA markers [4].The haplotype diversity of these two cpDNA fragments in the overall samples (the observed h = 0.9767 for the psbA-trnH fragment and the observed h = 0.8999 for the trnL-trnL fragment) is also greater than some other plants, such as Camellia huana (h = 0.759) [50], Michelia shiluensis (h = 0.674) [51], Iris loczyi (h = 0.820) [52], Rhododendron rex (h = 0.788) [53], Vouacapoua americana (h = 0.87) [54], and Tilia cordata (h = 0.881) [55].One main reason is the higher mutation rates occurring in these two fragments, although the general mutation rate over the whole chloroplast genome is not high [27].Many unique haplotypes with a couple of mutations that are different from major haplotypes likely occurred from recent mutations (Figure 1), which enhances the overall haplotype and nucleotide diversity.
The distribution of haplotypes in the western region was consistent with previous results.Xiao et al. [21] showed well-mixed phylogenetic relationships among individuals of four varieties of T. ciliata in Yunnan Province (T.ciliata var.ciliata, T. ciliata var.yunnanensis, T. ciliata var.pubescens, and T. ciliata var.henryi), implicating close genetic relationships among individuals in the western region in terms of chloroplast genomes.The distribution of haplotypes in the entire region showed an obvious phylogeographic structure, and different haplotypes were not randomly distributed in space.For example, dominant haplotype H1, coded in both fragments, appeared only in the eastern region of the whole distribution, and each was distributed in geographically similar ranges.These consolidate the previous study showing a phylogeographic structure in the T. ciliata complex with nuclear ITS markers [4].

Population Genetic Structure
The analysis of genetic structure among all populations showed moderate and significant genetic differentiation, accounting for 42.84% of the total genetic variation.This level of differentiation is lower than that derived from haploid nuclear SRAP (79.2%) [18], mtDNA (88.84%), and nuclear ITS (71.43%) [4] but greater than the genetic differentiation derived from SSR markers with the same 29 populations (32.4%) [20] and with six different populations (34.5%) [24].As expected, the genetic differentiation derived from cpDNA markers is greater than that of some tropical forest species derived from nuclear allozymes [56].In theory, when gene flow and genetic drift are the only two processes under the neutral hypothesis, population differentiation is expected to be greater for cpDNA markers of maternal inheritance than for nuclear markers of biparental inheritance but comparable between cpDNA and mtDNA markers with maternal inheritance [31].However, the preceding comparison is not in good agreement with this theoretical expectation, where genetic differentiation with cpDNA markers is smaller than that with ITS and mtDNA markers.One possible explanation is that high mutation rates in these two cpDNA fragments produced moderate genetic differentiation.This could be inferred in theory from Wright's formula under equilibrium of drift-migration-mutation in island model [30]: where m s is the migration rate of seeds and µ is the mutation rate.A high mutation rate results in small genetic differentiation, which was supported from a comparison of F st values derived from the SSR versus ITS markers (32.4% vs. 71.43%)with the same populations [4,20].The mutation rate of SSR markers is often considered to be higher than that of ITS markers.This could explain why the F st derived from cpDNA markers is smaller than that derived from mtDNA markers or potentially smaller than that derived from nuclear ITS markers despite lower genetic drift effects for nuclear genomes [31].A comparison of genetic differentiation between the two markers (Table 5) implies that the mutation rate in the psbA-trnH fragment could be greater than that in the trnL-trnL fragment, as was discussed by Zhu et al. [34] in medicinal plants.
A further insight into the effects of high mutation rates on genetic differentiation could be approximated from a comparison of F st(cp) = 1 + 2N e m s + µ cp −1 for cpDNA markers and F st(mt) = (1 + 2N e (m s + µ mt )) −1 for mtDNA markers under equilibrium.Assuming that the mutation rate for mtDNA markers is negligible (µ mt ≈ 0), we can obtain the following expression: Substitution of F st(mt) = 0.8884 [4] and F st(cp) = 0.4284 into Equation (3) yields a point estimate of µ cp /m s = 9.6234 ≫ 1.0, which is analogous to the genetic variation at copy number variation loci in Homo sapiens [57].The mutation rates at these two fragments are likely not as small as we generally thought before, and they play an important role in shaping the population genetic structure at these two loci.
Population structure in woody species is affected by multiple ecological and life history traits, including taxonomic status (angiosperm versus gymnosperm), regional distribution, geographical range, mating system, and seed dispersal [58].The geographical range of the T. ciliata complex involves complicated landscape structure (fragmented) and connectivity.From the results assayed with mtDNA markers [4], the investigated populations were divided into three groups from the elevations of 617.34 ± 264.19 m in the eastern region to 732.37 ± 250.87 m in the central region and 1407.67 ± 321.81 m in the western region.The central populations investigated (XY, LL, WM, XL, CH, TL, LD, and YF) were mainly located in Guizhou and Guangxi Provinces and showed small genetic variation among them due to geographic proximity.However, the overall range harbors diverse habitats that caused population isolation and intensified genetic differentiation among populations of the T. ciliata complex.Significant IBD effects could be an important factor causing substantial genetic differentiation at the range-wide scale, which was detected by both organelle and nuclear genome markers [4].The result indicating correlations between haplotype diversity and elevation implies that mountain ranges played a certain role in blocking gene flow at the range-wide scale [4].This is similar to the pattern of subtropical and tropical Machilus pauhoi in South China [59].The terrain in the eastern region is mostly mountain, which hinders both pollen and seed dispersal, resulting in less obvious pollen flow than that in the western region, in contrast to the western region, where pollen flow is more extensive than seed flow [4].
In addition, T. ciliata is mainly pollinated by insects [7], which could also enhance population structure.A study on the mating system of T. ciliata showed that pollen spread within and between populations was limited in different areas of the species' distribution [24].This may partly arise from the over-exploitation that led to low population density or the limited number of pollen donors within or between populations.Predominantly outcrossing with partial selfing and inbreeding facilitates IBD effects and, hence, population differentiation.The phylogenetic relationship among populations was closely linked to the geographical location, indicating an explicit haplotype geographical structure.

Figure 2 .
Figure 2. A map showing the twenty-nine sample sites and the geographic distribution of the cpDNA haplotypes of the psbA-trnH and trnL-trnL fragments: (A) the psbA-trnH fragment and (B) the trnL-trnL fragment.The pie charts show different proportions of haplotypes within each population.Different haplotypes are denoted by code Hi (i = 1, 2, … 47).Each color in the pie charts represents one haplotype (H1-H47).

Figure 2 .
Figure 2. A map showing the twenty-nine sample sites and the geographic distribution of the cpDNA haplotypes of the psbA-trnH and trnL-trnL fragments: (A) the psbA-trnH fragment and (B) the trnL-trnL fragment.The pie charts show different proportions of haplotypes within each population.Different haplotypes are denoted by code Hi (i = 1, 2, . . .47).Each color in the pie charts represents one haplotype (H1-H47).

Figure 3 .Figure 3 .
Figure 3. Tests of the effects of isolation by distance (IBD) on population genetic differentiation: (A) the psbA-trnH fragment and (B) the trnL-trnL fragment.Significant but weak IBD effects occurred Figure 3. Tests of the effects of isolation by distance (IBD) on population genetic differentiation: (A) the psbA-trnH fragment and (B) the trnL-trnL fragment.Significant but weak IBD effects occurred among populations.The red line in each figure indicated the trend of the relationship between Fst/(1-Fst) and the logarithm of geographic distance.

Figure 4 .
Figure 4. Phylogenetic relationships among 29 populations of the T. ciliata complex derived from concatenated sequences of the psbA-trnH and trnL-trnL fragments.The maximum likelihood method with 1000 bootstrap resamples was used to draw the consensus tree.Populations in grey and blue were grouped into eastern and western regions, respectively.

Figure 4 . 20 Figure 5 .Figure 5 .Figure 6 .
Figure 4. Phylogenetic relationships among 29 populations of the T. ciliata complex derived from concatenated sequences of the psbA-trnH and trnL-trnL fragments.The maximum likelihood method with 1000 bootstrap resamples was used to draw the consensus tree.Populations in grey and blue were grouped into eastern and western regions, respectively.Genes 2024, 15, x FOR PEER REVIEW 12 of 20

Figure 6 .
Figure 6.Population structure analysis of 447 individuals of the T. ciliata complex.(A) Cross-validation (CV) errors for the number of subpopulations (K) changing from 1 to 28. (B) A partitioned map of individuals with clustering assignments (K = 23) indicated in different colors.The proportions of genetic components from different subpopulations in each individual were indicated in different colors.

Figure 7 .
Figure 7. Analysis of mismatch distribution using the concatenated sequences of the psbA-trnH and trnL-trnL fragments: (A) population DC and (B) population TL.The grey bar lines and green lines are the observed frequencies under different numbers of segregation sites of pairwise sequences, whereas the red lines are the expected frequencies under a model of sudden population expansion [48].

Figure 7 .
Figure 7. Analysis of mismatch distribution using the concatenated sequences of the psbA-trnH and trnL-trnL fragments: (A) population DC and (B) population TL.The grey bar lines and green lines are the observed frequencies under different numbers of segregation sites of pairwise sequences, whereas the red lines are the expected frequencies under a model of sudden population expansion [48].

Table 1 .
• 14 ′ N, and 20.3 m above sea level) and the Yuejin North Experiment Station, South China Agricultural University (113 • 37 ′ E, 23 • 16 ′ N, and 42.3 m above sea level), in 2018.Note that some populations only had a few samples because supplementary samples were not obtained.All sample trees were separated by at least 50 m.The altitudes of the sampling sites ranged from 337 m (Guanshan, Jiangxi Province) to 1862 m (Huidong, Sichuan Province), with an average altitude of 808.2 m.Localities and 447 individuals of 29 populations of the T. ciliata complex analyzed with cpDNA markers.

Table 3 .
Haplotypes and haplotype diversity for the psbA-trnH fragment in each population of the T. ciliata complex.

Table 4 .
Haplotypes and haplotype diversity for the trnL-trnL fragment in each population of the T. ciliata complex.

Table 6 .
The neutrality test and mismatch distribution analysis with the concatenated alignment sequences of the psbA-trnH and trnL-trnL fragments in 29 populations of the T. ciliata complex.