The Genetic Diversity of Natural Ilex chinensis Sims (Aquifoliaceae) Populations as Revealed by SSR Markers

: Ilex chinensis Sims. is an evergreen tree species native to China and mainly distributed in the region south of the Qinling Mountains and the Huai River. This species has important ornamental, medicinal, ecological, and economic values, and plays a positive role in improving the environment and people’s lives. To reveal the genetic diversity and genetic structure of 401 individuals from 14 populations in the major distribution area of I. chinensis , 11 pairs of SSR primers were selected for PCR ampliﬁcation. The products were then subjected to capillary electrophoresis, and the ge-netic diversity of Ilex individuals was analyzed using relevant software. The results showed that the genetic diversity of I. chinensis was at a moderate-to-high level. A total of 54 alleles were detected at 11 SSR loci in the 14 Ilex populations, with an average of 4.831 alleles per locus. AMOVA analysis indicated that the genetic variation of I. chinensis populations mainly originated within populations. A STRUCTURE analysis divided the 401 I. chinensis individuals into four diﬀerent genetic clusters. The unweighted pair group methods using arithmetic averages (UPGMA) clustering based on Nei’s genetic distance revealed that the population from Xinping of Yuxi, Yunnan Province (XP), and the population from Longan of Qianxinan, Guizhou Province (LoA) were located in the outermost layer of the phylogenetic tree, indicating the furthest genetic relationship between these two population and other populations. The remaining populations could be roughly divided into two groups. Principal coordinate analysis (PCoA) demonstrated that the 401 individuals were clearly divided into three groups, which was consistent with the results of the STRUCTURE analysis and UPGMA clustering. This study identiﬁed the hotspots of genetic diversity of I. chinensis , as well as units for the conservation of individuals. It also revealed the pa�erns of genetic variation and population distribution of I. chinensis in diﬀerent regions, providing a molecular basis for the geographical zoning and formulation of breeding programs for I. chinensis , as well as germplasm resource management.


Introduction
Ilex, a genus in the monogeneric Aquifoliaceae, is unique phylogenetically and biogeographically.Due to its remarkable morphological diversity and phenotypic plasticity [1,2], this genus accommodates a vast number of species with inconsistent numbers by different authors, e.g., ≈400 spp.[3,4], ≈600 spp.[5], and >664 spp.[6].As a cosmopolitan genus, Ilex widely occurs in Asia, Oceania, Africa, Europe, the Canaries and Azores, North America, and South America, with an origin center in subtropical Asia [6] and modern diversification centers in East Asia and South America [7].
The role of Ilex for forest-community assembly and biodiversity maintenance is markedly significant in tropical and subtropical regions, as Ilex species are locally diverse in evergreen broadleaf forests.For instance, there are 14 and 12 species of hollies in the subtropical evergreen broadleaf forests of Mt.Wuyi and Mt.Huangshan in Eastern China, and four species (I.chinensis Sims., I. wilsonii Loes., I. editicostata Hu & Tang and I. suaveolens (H.Lév.) Loes.) are common in the la er [8,9].On Mt.Wuyi, I. elmerrilliana S. Y. Hu is a dominant species, with an important value of 2.64% [10].Ilex species may exhibit different life forms.In a mid-montane mixed forest on Mt.Emei in Western China, deciduous I. micrococca var.pilosa appears in the canopy layer, while evergreen I. fargesii Franch.and I. rockii S. Y. Hu appear in the small tree layer and the ground layer, respectively [11].For the Azores archipelago, with a temperate climate, I. azorica Gand. is commonly found dominant in natural forests [12].
In terms of utilization, holly seems to be versatile, with the most well-known being the source of beverages.In Chinese tea beverages, kudingcha (bi er spike tea) from holly leaves is very unique, and it is different from the tea from Camellia leaves.The leaf material of kudingcha was reported from several Ilex species, such as I. latifolia Thunb., I. kaushue S. Y. Hu (I.kudingcha), I. cornuta Lindl.& Paxton, and I. pentagona S. K. Chen, Y. X. Feng & C. F. Liang [13].I. paraguariensis A.St.-Hil.and I. guayusa Loes.from South America are sources of beverages due to their rich purine alkaloids (mostly caffeine, secondarily theobromine and theophylline) in leaves [14], and these beverages are called yerba maté and guayusa, respectively.Furthermore, hundreds of compounds have been identified from Chinese Ilex species, with the most common being triterpenoids, triterpenoid saponins, flavonoids, sterols, polyphenols, and other compounds [15].Thus, the Ilex species has multiple health benefits and potential for drug development [16,17].In addition, Ilex species are excellent ornamental and honey plants [18,19].
Great a ention has been paid to the species richness, geographical diversity, ecological importance, and unique utilization of Ilex species.Therefore, a number of biological investigations have been focused on taxonomic, phylogenetic, and genomic studies.On taxonomy, new taxa continue to be discovered [20][21][22].In terms of phylogeny, DNA molecular evidence has been widely used to reveal evolutionary relationships between species.A series of regional and global studies have been reported [1,[4][5][6][7]23,24].As for genome studies, currently, only a small number of reports focus on the whole genome [25,26] and mitochondria genome [27].More reports are focused on chloroplast genomes [28][29][30].Intraspecific microevolution, however, has not received sufficient a ention yet.The reports on genetic diversity and genetic structure are limited to a few species, such as I. paraguariensis [12,[31][32][33][34], I. guayusa [35], and I. crenata Thunb.[36].Although molecular markers have developed for I. leucoclada (Maximowicz) Makino [37] and I. chinensis [38], there is no evaluation of genetic variation.For some species, genetic diversity has been reported; however, the markers used are not codominant [39].In addition, in Europe, there is a report on a chloroplast haplotype analysis of the English holly species (I.aquifolium Marshall) [40].
In summary, investigations into the genetic variation of Ilex tend to focus on the super-specific level (mostly based on small samples for a species) and neglect the intra-specific level (based on populations).For a given holly species, an in-depth understanding of the pa erns of genetic variation among populations and geographical regions, as well as the driving factors that contribute to it, is still lacking.This situation is not informative to the establishment of the species concept for hollies.Because a phenomenon of hybridization and introgression seems to occur in some hollies [41][42][43], this results in blurred species boundaries.On the other hand, insufficient information on genetic variation at the population level limited the practice in the collection, preservation, and utilization of Ilex species germplasm resources.
As an evergreen tree species, I. chinensis is native to China and mainly occurs in the area south of the Qinling Mountains and Huai River, usually found in evergreen broadleaved forests and forest edges on mountain slopes at an altitude of 500-1000 m.In addition to China, it has also been extensively cultivated and introduced to other regions around the world, such as Japan, South Korea, Vietnam, Malaysia, and others.This species has significant ornamental, medicinal, ecological, and economic value, and plays a positive role in environmental improvement and people's lives.There is only one report on the development of molecular markers for genetic variation for I. chinensis, and it lacks a systematic analysis of genetic diversity.As a widely distributed and evergreen tree species, it is an ideal material for the exploration of genetic variation of Ilex species.The aim of this study is to (i) explore natural populations of I. chinensis, (ii) analyze intraspecific genetic diversity using microsatellite markers, (iii) reveal the distribution pa erns of genetic variation between populations and regions, and (iv) infer potential environmental factors contributing to genetic variation pa erns.

Population Sampling
From August 2022 to April 2023, 401 individuals were sampled from 14 natural populations of I. chinensis in China (Figure S1).In each population, 10-15 fresh leaves per tree were collected from 20-30 adult individuals.In order to avoid collecting individuals with the same genetic background, the distance between adjacent individuals in the same population should not be less than 30 m.For smaller populations, all individuals were collected.Leaf tissues were quickly dried using silica gel and stored at room temperature in the laboratory.The geographical locations and altitudes of the sampled populations were recorded using the Biotracks app (Table S1).

DNA Extraction and Microsatellite Genotyping
DNA extraction was performed using the Rapid Plant Genome DNA Extraction System (DP321) Kit (Tiangen, Beijing, China).The concentration and purity of the extracted DNA were measured by a NanoReady Micro UV-Vis spectrophotometer (Life Real, Hangzhou, China), and the DNA was diluted to an appropriate concentration with TE buffer.For the qualification of the DNA sample, 0.8% agarose gel electrophoresis was conducted, and the qualified samples were stored in the refrigerator at −20 °C.
Based on the transcriptome sequencing results of Ilex × altaclerensis 'Belgica Aurea' [44], 80 SSR loci primers were picked and synthesized by the General Biol company (Chuzhou, China), and 11 SSR loci primers suitable for I. chinensis were selected from the 80 SSR primers to study the genetic diversity of the I. chinensis genotypes (Table S2).
The 5′ ends of the forward sequences of the screened polymorphic primers were fluorescently modified (FAM, HEX, TAMRA, or ROX).All samples were amplified using fluorescence-modified forward primers and normal reverse primers, and polymerase chain reactions (PCRs) were performed with 2×Taq PCR MasterMix (Tiangen, Beijing, China) with an amplification system of 20 µL.The PCR protocol was designed as 35 cycles of initial denaturation at 94 °C for 4 min, followed by an extension of 35 s at 94 °C, appropriate annealing temperature for 35 s, extension at 72 °C for 45 s, and ending with a final extension at 72 °C for 5 min.PCR amplification products of different lengths and different fluorescent markers were mixed for genotyping by a high-throughput genetic analyzer (ABI3730, Applied Biosystems, Waltham, MA, USA).

Data Analysis
Sequencing results from capillary electrophoresis were imported into GeneMarker 2.2 software (SoftGenetics, State College, PA, USA), and parameters such as length intervals of molecular internal standards and target fragments were set to read the band length data of microsatellite primer amplification products.
The number of observed alleles (Na), the number of effective alleles (Ne), the Shannon information index (I), the observed heterozygosity (Ho), the expected heterozygosity (He), the genetic differentiation index (Fst), gene flow (Nm), analysis of molecular variance (AMOVA), and principal coordinate analysis (PCoA) were determined using GenAlEx 6.503 software [45].The polymorphism information content (PIC) was determined using PowerMarker software (version V3.25) [46].Based on Nei's standard genetic distance, MEGA7 software was used to perform a systematic clustering analysis on the populations [47].A Bayesian clustering analysis in STRUCTURE 2.3.4 software was used to analyze populations' genetic structure and estimate individual ancestry [48].The map was drawn using ArcGIS 10.8.2 software [49].

Genetic Diversity at the Species Level
A total of 54 alleles were detected at 11 SSR loci in 14 I. chinensis populations, with an average of 4.831 alleles per SSR locus (Figure 1).Na per locus ranged from 2.000 (Ilex-SSR-32) to 10.786 (Ilex-SSR-75).Ne was an indicator reflecting the size of genetic variation within a population.The closer this value was to the detected Na, the more evenly distributed the alleles were in the population.The mean Ne in this study was 2.573, which significantly differed from the mean Na.It suggested that the distribution of alleles in I. chinensis populations overall was not uniform enough.Ho ranged from 0.037 (Ilex-SSR-32) to 0.796 (Ilex-SSR-75), with a mean of 0.440.He ranged from 0.09 (Ilex-SSR-32) to 0.822 (Ilex-SSR-75), with a mean of 0.488.In plant-based studies, when He was around 0.61 ± 0.21, it indicated that the genetic diversity of the species was at a moderate to high level.Therefore, it could be seen that the I. chinensis population had a relatively high genetic diversity.Similarly, Shannon's information index (I) ranged from 0.181 (Ilex-SSR-32) to 1.991 (Ilex-SSR-75), with a mean of 0.977, indicating a relatively high level of genetic diversity in I. chinensis at the species level.
In this study, all values of F obtained were positive, suggesting that homozygotes were dominant in the population.For I. chinensis, Fis ranged from 0.008 (Ilex-SSR-01) to 0.590 (Ilex-SSR-32), with a mean of 0.163, indicating the possibility of inbreeding within local populations, but at a relatively low frequency.On the other hand, Fit ranged from 0.123 (Ilex-SSR-75) to 0.795 (Ilex-SSR-32), with a mean of 0.333, indicating a relatively low overall frequency of inbreeding within the entire population.These results were consistent with the characteristics of I. chinensis as an entomophilous plant with long pollendispersal distances, weak self-compatibility, and predominantly single-plant distribution, which may have resulted in a relatively distant genetic relationship between parents.
According to Table 1, it can be seen that among the 11 microsatellite loci selected, only Ilex-SSR-32, Ilex-SSR-34, and Ilex-SSR-41 exhibited moderate-to-low polymorphism, while all other loci were highly polymorphic.The average PIC for each locus was 0.587, ranging from 0.176 (Ilex-SSR-32) to 0.911 (Ilex-SSR-75).This indicated that the microsatellite loci selected in this study had an overall high polymorphism.

Genetic Diversity at the Population Level
The genetic diversity indices for 14 natural populations of I. chinensis are shown in Table 2.The range of Na was between 3.636 and 6.182, with an average of 4.831.Ne ranged from 2.198 to 2.939, with an average of 2.573.The range of Shannon's information index was between 0.829 and 1.198, with an average of 0.927.Ho ranged from 0.365 to 0.535, with an average of 0.440.The minimum value was for the WS population, and the maximum was for the LiA population.He ranged from 0.402 to 0.595, with an average of 0.488.The minimum value was for the JS population, and the maximum was for the JR population.Among the 14 populations of I. chinensis, the He values of five populations were greater than 0.5, indicating that these populations had relatively rich genetic diversity.The Ho value of the LY population was greater than the He, indicating a low level of homozygote excess in this population.In the remaining 13 populations, Ho values were all smaller than He, indicating an excess of heterozygotes in these populations.Based on these comprehensive genetic diversity indices, the JR population had the highest genetic diversity, while the JS population had the lowest genetic diversity.
From a geographical perspective, there were some differences in the genetic diversity of the I. chinensis populations in the southeastern region, central region, and southwestern region of China.In particular, He for the southeastern region was 0.515, the central region was 0.483, and the western region was 0.461.Therefore, the overall genetic diversity of I. chinensis in the southeastern region of China was slightly higher than that in the central and southwestern regions of China.

Genetic Variation and Genetic Structure
The F-statistic (Fst) was 0.185 (Table 3), which fell within the range of 0.15 to 0.25.It indicated that there was a high level of genetic differentiation among the I. chinensis populations, suggesting that there were substantial differences between them.This might be due to the relatively low average gene flow between populations (Nm = 0.870, Table 1).An AMOVA analysis of 14 different I. chinensis populations indicated that 18.5% of the genetic variation was found among populations, while 81.5% was found within populations (Table 3).This suggested that there was a high level of genetic variation within individual populations and that the main source of genetic variation in I. chinensis came from within-population processes.
In order to analyze the genetic structure of I. chinensis populations, a Bayesian clustering model was used to analyze the ancestral relationships among these populations.The results of the ΔK analysis showed that the maximum value of ΔK was reached when K = 4 (Figure 2).The clustering result for K = 4 is shown in Figure 3, with four different colors (green, blue, yellow, and red) representing four genetic clusters named qI~qIV (Table S3), respectively.According to the geographical distribution map of the genetic structure of I. chinensis populations (Figure 4), most of the genetic information of populations WZ, LY, JR, and LiA in the southeast came from the same ancestral population.Similarly, populations located in the central region and some parts of the southwest also shared most of their genetic information from the same ancestral population.The populations JZ, LiA, and WYS, located in the Qinling Mountains, Dabie Mountains, and Wuyi Mountains, respectively, showed obvious gene introgression.In contrast, the XP and LoA populations in the southwest may have had less gene flow due to their higher altitudes, which led to a greater genetic distance from other populations.A UPGMA tree of 14 I. chinensis populations was constructed based on Nei's genetic distance (Figure 5), in which XP and LoA were located at the outermost part of the phylogenetic tree, indicating that this population had the farthest genetic relationship with other populations.The remaining population could be roughly divided into two groups at a genetic distance of 0.1.Among them, the four southeastern populations of JR, LiA, WZ, and LY clustered into one branch, while YL, PiJ, SC, WS, JZ, WYS, JS, and PuJ clustered into another branch.This branch included all central populations and some southeastern and southwestern populations.The clustering analysis results indicated that there was a clear geographical structure between populations, meaning that most populations with closer geographical locations gathered together first, which was similar to the results of STRUCTURE.A principal coordinate analysis (PCoA) based on the inter-population pairing Fst matrix of the 14 I. chinensis populations revealed that the 401 individuals were clearly divided into three groups (Figure 6).Group 1 included 12 populations distributed in the southeastern, central, and some southwestern regions.Group 2 consisted of the XP population located in the southwestern region, while Group 3 included the LoA population in the southwestern region.These results were consistent with the findings of the structure analysis and similar to the results of the UPGMA clustering analysis.

Genetic Diversity
The SSR molecular markers are crucial for genetic diversity analysis, due to their several genetic features, such as codominance inheritance [50].SSR markers may have higher polymorphism than dominant markers (RAPD, ISSR, etc.).As a key index of a plant's adaptability to environmental changes, genetic diversity may be influenced by natural selection, genetic drift, migration, and reproductive systems [51].Studying the genetic diversity of Ilex at the population level will provide a theoretical reference for the protection and utilization of Ilex species resources.Although the species in Ilex are important types of resource plants, with high ecological, ornamental, and economic value, research on the genetic diversity of this group is currently weak, especially for I. chinensis.
Here, for the first time, we provide an empirical study of the genetic diversity across 14 populations of I. chinensis.There are not many comparable case studies of genetic diversity for Chinese hollies; AFLP Markers were used to analyze the genetic diversity of one wild and five cultivated populations of I. glabra (L.) A. Gray for the identification of the varieties, genetic improvement, and germplasm resource protection [52].The development of SSR markers for I. chinensis has been reported [38], but no subsequent diversity analysis was found.ISSR markers have been used for I. integra Thunb., an insular endemic species of hollies, indicating a low genetic diversity and a high genetic differentiation [39].It is obvious that the measured data of SSR markers for Chinese hollies is basically lacking.As for the Ilex species materials from other regions, there are several genetic diversity studies available for reference.
Investigation of six populations from successive age cohorts (adult, regenerating, and seed-derived populations) of I. paraguariensis in Southern Brazil revealed a relatively high level of He, with the highest and lowest values of 0.660 and 0.584, respectively [34].The overall He across 14 nuclear SSR loci was 0.459 for the "yerba mate" germplasm of I. paraguariensis [31].Another report also showed that "yerba mate" breeding populations were highly variable [33].On the other hand, the results from I. guayusa suggested a moderately low degree of genetic diversity (He = 0.396) [35].
Based on our data, the overall He of the 11 SSR loci for I. chinensis was 0.488, indicating that this holly species had a moderately high genetic diversity.There are slight differences in He among different geographical groups, with an average He of 0.515, 0.483, and 0.461 for the southeast (five populations), central (five populations), and southwest groups (four populations), respectively.The trend of He between geographical groups is also supported by other genetic parameters, and there is a downward trend from the east to west geographical groups for the number of efficient alleles (Ne), the observed heterozygosity (Ho), and Shannon's information index (I).Overall, compared to existing reports, I. Chinensis holds a moderately high level of genetic diversity and shows a decreasing trend in geographical groups from east to west.
UPGMA clustering, STRUCTURE clustering based on the Bayesian algorithm, and PCoA all had similar results, which means that those with similar geographic locations were clustered together.Regional genetic clustering was based on population, specifically origin clustering, and this genetic structure was typical in isolated populations.Among them, the XP and LoA populations located in high-altitude areas constitute two separate branches.All southeast populations, except population WYS, form a genetic cluster with a central cluster assembled by other populations.The geographical pa ern of this genetic structure roughly follows the overall trend of east-central-west differentiation, although there are exceptions.In fact, previous reports are generally consistent.Castanopsis fargesii Franch.was reported to display a typical geographic pa ern of east-central-west [61].An obvious pa ern of the two clades was reported for the lacquer tree, e.g., a western group including 35 sites in Western China, and the eastern group embracing 40 sites in Eastern China, South Korea, and Japan [62].Three genetic clades (eastern, southwestern, and northwestern regions) were found in Liriodendron chinense (Hemsl.)Sarg.[63].
Overall, I. Chinese has a relatively high genetic differentiation among populations, and the geographical pa ern of genetic structure is manifested as three clades: east, central, and west.

Intrinsic Factors Affecting Genetic Variation
Life histories and breeding systems highly affect genetic diversity and structure [64].As a group of typical dioecious species, Ilex species have small and white or pale-colored flowers and are pollinated by insects [65].So far, there have been no reports on the breeding system and pollen and seed vectors of I. chinensis.Nevertheless, based on its dioecious sex system and relevant observations of the hollies mentioned above, it is not difficult to speculate that the breeding system of I. chinensis is outcrossing.Apis cerana Fabricius was the flower visitor of fourteen holy species [66].A recent investigation has shown that βcaryophyllene in I. rotunda Thunb.may be used as an information substance to a ract Apis cerana for pollination [67].In a dioecious species, I. opaca Aiton, female reproductive success and fruit yields were limited by resource levels rather than pollen availability [68].Due to its dioecy, holly species typically maintain a breeding system of outcrossing.Outcrossing species tended to be more genetically diverse and had less genetic differentiation among their populations, while woody plants had less population differentiation [64].In addition, the mode of reproduction (sexual vs. asexual) is likely to have important effects on genetic variation and its spatial distribution within plant populations [69,70].Based on relevant expectations [64], outcrossing can explain why I. chinensis, to a considerable extent, has moderately high genetic diversity, rather than relatively high genetic differentiation.
Hollies were reported to be bird-dispersed tree or shrub species, such as I. integra [71] and I. aquifolium [72,73].Although the main mode of seed transmission is bird dispersal, such as Turdus iliacus L., the seeds are secondary transported by predating rodents, such as Clethrionomys gladiolus [74].Interestingly, the seeds of I. aquifolium may also be spread by mammals (mainly Vulpes vulpes and Meles meles) [75].The phenomenon of post-dispersal seed predation also exists in I. canariensis Poir.[76].In another study, different fruits of I. chapaensis Merr.were documented by various animals, with Paguma larvata for large green fruits and birds for red or black fruits [66].Seed dispersal by birds is assumed to profoundly impact genetic diversity in many plant species [77].Environmental fecal samples in a study of plant-frugivore interactions can be applied to the investigation of plant ecology and genetic variability [78].Although long-distance dispersal syndromes (ornithochory) do not always ensure genetic connectivity across large areas in island systems [79], we believed that it was crucial to promote gene flow among populations of I. chinensis.Thus, relatively moderate genetic differentiation between populations may be related to the obstacle of plant seeds-birds interaction.

Environmental Factors Affecting Genetic Variation
As for forest plants, altitude has a potential impact on genetic diversity within species.However, a previous inquiry demonstrated different pa erns of variation.Populations at intermediate altitudes have greater diversity than populations at lower and higher altitudes [80].The data of genetic diversity analysis from five subpopulations of Rhododendron simsii Planch.at altitude gradients (972-1370 m) showed a high-low-high variation pa ern in expected heterozygosity (He) [62].Contrasting elevational pa erns of genetic variation were reported in Euptelea pleiospermum Hook.f. et Thoms., and significant impacts were found for the central mountains, with higher genetic diversity in the middlealtitude populations than in the low and high altitudes.No elevational effects were detected for marginal mountains [81].In a study on grassland species, the authors drew a conclusion that altitude does not affect genetic diversity [82].For the case of I. chinensis, the altitude variation of genetic diversity presented a negative correlation, with high diversity in low-altitude populations and low diversity in high-altitude populations.Specifically, the populations in the southeast were widely distributed in the Middle and Lower Yang e Valley Plain, while the populations in the southwest were only sca ered in the high-altitude areas of the Wumeng Mountains and Hengduan Mountains.
Ilex was believed to have originated in the early Eocene (50.8 Ma), grown in moist subtropical East Asia, and expanded its range of temperature tolerance rather than moisture tolerance in the subsequent evolutionary process [6].I. chinensis was dated from the Pliocene-Pleistocene boundary (2.53 Ma) [6].The uplift of the Himalaya-Hengduan Mountains (late Miocene to early Pliocene) accompanied by simultaneous cooling was assumed to trigger the differentiation of oaks [83].Specifically, the uplift of the Yunnan-Guizhou Plateau began at 2.8 Ma and had a fast and dramatic uplift around 0.8 Ma [84].Thus, we speculated that the eastern-central-western pa ern of genetic differentiation for I. chinensis may be affected by the uplift of the Yunnan-Guizhou Plateau and its associated climate change.In addition, human interference may also affect the genetic structure of I. chinensis.The fine-scale genetic structure of I. aquifolium may be influenced by human interventions (mainly coppice practice and ca le browsing) [85].
Population JR appeared to be mixed with a few individuals similar to the individuals of population XP (Figure 3).We do not rule out the existence of introduced individuals in population JR.However, further exploration is needed.

Conclusions
This study found for the first time that the genetic diversity of I. chinensis was at a moderate level, and most of this genetic variation existed within populations.The high level of population differentiation could have been maintained through bird-mediated seed dispersal.The geographic pa ern of I. chinensis showed a gradual decrease in genetic diversity from east to west and from low to high altitudes.This study provided the complete genetic data of I. chinensis for the first time, which provided a molecular basis for the collection of genetic resources, the preservation of core germplasm resources, and the development of new varieties of I. chinensis.It also provided valuable information for hybridization and species boundary delineation of I. chinensis.For deeper genetic diversity analysis, a series of follow-up studies, such as cpDNA and nDNA sequencing, genomics, etc., would be needed.

Supplementary Materials:
The following supporting information can be downloaded at h ps://www.mdpi.com/article/10.3390/f15050763/s1: Figure S1: The geographical distribution of 14 populations of Ilex chinensis; Table S1: Accession data for the populations of Ilex chinensis; Table S2: Primer sequences and fragment size information for 11 microsatellite loci; Table S3: The proportion of four genetic branches in 14 populations.
Author Contributions: Conceptualization, Q.Z. and M.Z.; formal analysis, X.W.; investigation, S.H. and P.Z.; resources, S.H. and P.Z.; data curation, X.W.; writing-original draft preparation, S.H. and Y.F.; writing-review and editing, Q.Z. and M.Z.; funding acquisition, Q.Z. and M.Z.All authors have read and agreed to the published version of the manuscript.

Figure 1 .
Figure 1.Maps of PCR product separation for 11 SSR loci.The first three individuals (A-C) came from the WZ population, and the last three individuals (D-F) came from the SC population.

Figure 2 .
Figure 2. Relationships between the number of clusters (K) and the corresponding ΔK statistics calculated according to ΔK (Evanno et al., 2005) [48] based on Structure analysis.

Figure 3 .
Figure 3. Population genetic structure based on the Bayesian clustering model among 216 samples at K = 4. Green, blue, yellow and red represent the four genetic clusters of qI~qIV, respectively.

Figure 4 .
Figure 4. Geographical distribution of STRUCTURE results and genetic structure in 14 populations (K = 4).The WZ, LY, JR, LiA, and WYS populations belong to the southeastern population.The JZ, JS, YL, PiJ, and SC populations belong to the central population.The WS, PuJ, XP, and LoA populations are from the southwest.Green, blue, yellow and red represent the four genetic clusters of qI~qIV, respectively.

Figure 6 .
Figure 6.Score plot generated using PCoA.All the results showed that the samples could be divided into three clusters.
Na: number of alleles, Ne: effective number of alleles, I: Shannon information index, Ho: observed heterozygosity, He: expected heterozygosity, F: fixation index, Fis: inbreeding coefficient of a local population, Fit: overall inbreeding coefficient, Fst: F-statistic, Nm: geneflow, PIC: polymorphism information content.

Table 2 .
Mean values of the genetic diversity statistics for 11 microsatellite loci in 14 populations.