Expressed sequence tag-simple sequence repeats were identified in transcript sequences from the Sapindus whole transcriptome sequencing project using MISA¹ (Thiel et al., 2003), and Primer3 (Untergasser et al., 2012) was used to design EST-SSR primers. We screened all EST-SSR markers against eight Sapindus germplasm accessions representing different geographic origins and selected 18 pairs of highly polymorphic and stable EST-SSR markers (Table 1) (unpublished). The 18 EST-SSR primers were synthesized by Beijing Ruiboxingke Biotechnology Co. Ltd. and used in subsequent experiments.

The Sapindus population analyzed in this study comprised 161 wild individuals, with 160 individuals from 16 provinces in China and one individual from Vietnam (Figure 1 and Supplementary Table 1). It included 117 S. mukorossi individuals, as well as 36 S. delavayi, four S. rarak, and four S. rarak var. velutinus representatives. The criteria for Sapindus germplasm collection were representative local, naturally superior plants with a diameter at breast height of at least 30 cm. These germplasms were conserved ex situ by grafting at a Sapindus national germplasm nursery in Jianning County, Fujian Province, China (27°06′ N, 117°25′ E), comprising 175 clones and 64 half-sib families by grafting and seeding respectively, with at least 15 plants for each clone. The average temperature in germplasm nursery is 17.4°C, with a maximum temperature of 36.8°C (July) and a minimum temperature of 4.3°C (January), and an average relative humidity of 83.9% (Wang et al., 2020). At present, the grafted clones are 6 years old.

300 mg fresh leaf tissue from each Sapindus accession was used for DNA extraction with a Plant Genomic DNA Extraction Kit (Tiangen, DP320-03). The DNA samples concentration and quality were determined with a VERTEX 70 UV spectrophotometer (Bruker, Germany).

Eighteen EST-SSR markers were applied to genotype the 161 Sapindus tree samples. The forward primer of each marker was labeled at the 5′ end with a fluorescent dye (FAM, HEX, TAMRA, or ROX) according to the target fragment size of the marker. PCR (Polymerase Chain Reaction) was performed on a T100 thermal cycler (Biorad) in a 20-μL volume that contained 10.0 μL 2×Taq Plus PCR mix [Taq Plus DNA Polymerase (recombinant), 3 mM MgCl₂, 0.4 mM dNTPs (dATP, dCTP, dGTP, and dTTP): 0.4 mM], 0.3 μL 10 μM forward primer, 0.3 μL 10 μM reverse primer, and 1.0 μL genomic DNA. The PCR protocol was 5 min denaturation at 94°C; 20 cycles of 30 s at 94°C, 30 s at 52/72/95/50°C (annealing temperature depended on the primer used, see Table 1), and 30 s at 72°C; and a final extension for 5 min at 72°C. Microsatellite alleles were called on an ABI 3730XL DNA analyzer, and the amplicons were statistically analyzed by GeneMarker 2.2.0 software (SoftGenetics, State College, PA, United States).

For each microsatellite marker, the number of alleles (Na), number of effective alleles (Ne), observed heterozygosity (Ho), expected heterozygosity (He), unbiased expected heterozygosity (uHe), ibreeding coefficient (F), and Shannon’s information index (I) were calculated using GenAlEx v6.5 (Peakall and Smouse, 2006). The polymorphic information content (PIC) was calculated using PowerMarker V3.25 software (Liu and Muse, 2005). F-statistics calculations (FIS, FIT, and FST) and principal coordinate analysis (PCoA) were also performed in GenAlEx v6.5 combined with Microsoft Excel. A neighbor-joining (NJ) tree was generated based on pairwise genetic distances between individuals by using PowerMarker V3.25 (Liu and Muse, 2005), and plotted with iTOL.²

The population structure analysis was performed with Bayesian model-based admixture analyses in STRUCTURE v2.3.4 (Pritchard et al., 2000). We set the number of Markov chain Monte Carlo (MCMC) iterations after burn-in to 100,000 with a 100,000-run length, and the number of genetically homogeneous clusters (K value) ranged from 1 to 20 with 10 replicate runs for each analysis. The optimum K-value was determined by the highest ΔK method (Evanno et al., 2005) in Structure Harvester³ (Earl and VonHoldt, 2012). The structure plot was constructed in R 4.1.0 (R Core Team, 2013).

The optimal set of core germplasm was extracted by the Core Hunter 3 (De Beukelaer et al., 2018) which maximized the genetic variation and allelic richness using local search algorithms. Based on the previously reported distribution of core germplasm fractions in woody plants ranging from 10 to 45% (Belaj et al., 2012; Duan et al., 2017; Feng et al., 2018; Lv et al., 2020), we decided to test 10 sampling fractions (10, 15, 20, 25, 30, 35, 40, 45, and 50% and initial group) respectively by Core Hunter 3. Na, Ne, Ho, He, I, and uHe were calculated separately for each fraction using GenAlEx software, as described above. These indicators were t-tested between the core subset and the initial group using Microsoft Excel. The smallest core subset that did not differ significantly with the 100% population group (P ≤ 0.05) was then selected as the optimal core germplasm collection (Lv et al., 2020).

There were 236 alleles identified by the 18 EST-SSR markers. All 161 individuals could be uniquely genotyped using these 236 alleles, demonstrating the high discrimination capacity of these 18 EST-SSR markers. The markers showed considerable variation (Table 2), with number of alleles (Na) ranging from 6 to 25 (mean 13.1), number of effective alleles (Ne) from 2.625 to 11.503 (mean 5.711), observed heterozygosity (Ho) from 0.277 to 1.000 (mean 0.558), expected heterozygosity (He) from 0.619 to 0.913 (mean 0.798), unbiased expected heterozygosity (uHe) from 0.621 to 0.916 (mean 0.801), Shannon’s information index (I) from 1.202 to 2.683 (mean 1.937), and polymorphic information content (PIC) from 0.561 to 0.907 (mean 0.775). The marker with the highest number of alleles was S36 (25), and the marker with the lowest was S140 (6). The observed heterozygosity of all markers was lower than the expected heterozygosity, with the exception of S105. All markers exhibited high polymorphism (PIC > 0.5; Table 2). The S36 marker captured the most genetic diversity with the highest PIC value (0.907).

There were marked differences in genetic diversity among the four different Sapindus taxa (Table 3). S. mukorossi exhibited the highest Na, Ne, I, Ho, He, and uHe values and the lowest F value compared with S. delavayi, S. rarak, and S. rarak var. velutinus. By contrast, S. rarak var. velutinus showed lower levels of genetic diversity, probably because of its narrow distribution area and smaller number of individuals. The mean pairwise F_ST coefficient between the four species of Sapindus was 0.154 (Table 4); highest differentiation was between S. delavayi and S. rarak (0.183) and lowest for S. mukorossi versus S. delavayi (0.122).

The first and second coordinates of the PCoA analysis accounted for 40.54 and 13.44% of the total genetic variation, respectively (Figure 2). The population of Sapindus was split into three clusters on the first principal coordinate axis (PCoA 1), which corresponded approximately to the S. delavayi group, the S. rarak and S. rarak var. velutinus group, and the S. mukorossi group. On the second principal coordinate axis (PCoA 2), individuals of S. mukorossi were roughly divided into two subgroups: the individuals from southern Guizhou province and the remaining S. mukorossi germplasm. Likewise, S. rarak and S. rarak var. velutinus were also divided into two subgroups along the second principal coordinate axis.

In contrast to the PCoA results, the STRUCTURE analysis indicated that the Sapindus population could be genetically divided into two distinct subgroups by the ΔK method (Figures 3B,C). Subgroup 1 contained all individuals of S. mukorossi, and subgroup 2 included all individuals of S. delavayi, S. rarak, and S. rarak var. velutinus. When K was equal to 3–6, individuals of subgroup 2 were consistently divided into S. rarak subgroup and S. rarak var. velutinus subgroup (Figure 3A). Furthermore, S. mukorossi individuals were divided into several subgroups (K = 3–6). Intriguingly, individuals from southern Guizhou province formed a separate subgroup (K = 3–6), suggesting that they originated from a distinct ancestral population.

The neighbor-joining dendrogram based on genetic distances among individuals (Figure 4) grouped S. mukorossi, S. delavayi, and S. rarak into their own populations, except for one S. delavayi individual that was grouped together with the S. mukorossi individuals. Furthermore, S. mukorossi individuals could be divided into multiple subgroups which largely correlated with their geographic distribution. Interestingly, some S. mukorossi individuals from Guizhou province appeared to be more closely related to S. rarak.

We produced nine candidate core germplasm collections of different sizes using Core Hunter. To determine the optimal core germplasm size, we compared the nine subgroups to the whole population according to six diversity parameters: Na, Ne, Ho, He, uHe, and I (Table 5). The core 10 and 15% subgroups differed significantly (P ≤ 0.05) and highly significantly (P ≤ 0.01) from the full population in Ne and I, respectively. The core 10–25% subgroups differed highly significantly (P ≤ 0.01) and significantly (P ≤ 0.05) in Na from the core 30 and 35% subgroups. The remaining parameters showed no significant differences between the subgroups and the full collection. Hence, the core 40% subgroup was selected as the optimal core germplasm collection. It contained 64 Sapindus individuals: 53 S. mukorossi, nine S. delavayi, one S. rarak, and one S. rarak var. velutinus (Supplementary Table 2), respectively. Relative to the full collection, the core 40% subgroup maintained 80.66% of Na, 97.65% of I, and 90.69% of Ho, and it had 101.90% higher Ne, 103.68% higher He, and 109.59% higher uHe.

The collection of natural germplasm resources and genetic diversity evaluation is important for the conservation, breeding, and utilization of germplasm resources (Glaszmann et al., 2010). Previous studies have documented abundant genetic diversity in Sapindus. Diao et al. (2016) and Jiang et al. (2016) found significant genetic divergence among germplasm accessions of S. mukorossi using ISSR markers. Sun et al. (2018) also found significant interspecific genetic differences between S. mukorossi and S. delavayi in China using ISSR markers, and they reported that S. mukorossi could be broadly divided into two subgroups. In this study, for the first time, we have pooled S. mukorossi, S. delavayi, S. rarak, and S. rarak var. velutinus germplasms from China to analyze their genetic diversity and population structure using EST-SSR markers. We also found substantial genetic diversity within the Sapindus germplasm, with S. mukorossi exhibiting the highest genetic variation (Na = 11.83, I = 1.82); S. rarak and S. rarak var. velutinus showed less variation (Na = 2.94 and 2.89 and I = 0.83 and 0.85, respectively) (Tables 2, 3). This may be due to the wider distribution and population size of S. mukorossi; S. rarak and S. rarak var. velutinus have a limited distribution in Yunnan Province, China, with S. rarak found only in Xishuangbanna Dai Autonomous Prefecture, Yunnan Province and S. rarak var. velutinus only in Shiping and Jianshui counties in Honghe Hani and Yi Autonomous Prefecture, Yunnan Province.

It is difficult to distinguish S. mukorossi, S. delavayi, S. rarak, and S. rarak var. velutinus in nature because of their high phenotypic similarity, and there has been a lack of molecular biological support for the species splits in the Sapindus taxonomy. A better understanding of Sapindus genetic structure is a first step toward addressing these issues. Here, we have identified varying levels of genetic divergence among S. mukorossi, S. delavayi, S. rarak, and S. rarak var. velutinus for the first time using SSR markers. Pairwise F_ST comparisons among the four Sapindus taxa showed high genetic differentiation (Table 4), and PCoA (Figure 2) and a neighbor-joining dendrogram (Figure 4) supported this (Figure 2). However, structure analysis divided the Sapindus germplasms into two subgroups, a S. mukorossi subgroup and subgroup containing the other taxa (Figure 3A), probably due to the higher representation of S. mukorossi among the samples. Structure analysis was developed for comparing populations within a single species and requires the assumption of Hardy-Weinberg equilibrium in the population, which may be violated by the populations of Sapindus studied here (Lv et al., 2020). Hence, PCoA analysis can be a more valid and efficient approach for genetic structure identification owing to its relaxed Hardy-Weinberg equilibrium hypothesis (Lv et al., 2020). We presented the results of structure analysis for K = 3–6 and found that S. mukorossi, S. delavayi, and the combination of S. rarak and S. rarak var. velutinus were divided into three distinct subgroups; S. mukorossi could be further divided into multiple subgroups (Figure 3A). Intriguingly, structure analysis suggested that S. rarak and S. rarak var. velutinus originated from the same ancestral population (Figure 3A), and PCoA analysis also indicated that S. rarak and S. rarak var. velutinus were closely related (Figure 2), supporting the theory that S. rarak var. velutinus is a variety of S. rarak. Moreover, we also found that individuals from southern Guizhou Province formed a distinct subgroup (Figure 3A), consistent with the PCoA analysis in which S. mukorossi individuals were divided into two subgroups along PCoA axis 2 (Figure 2). Previously, (Mahar et al., 2011a,b, 2013) used RAPD, DAMD, and ISSR molecular markers to analyze germplasm of S. mukorossi, Sapindus trifoliatus, and Sapindus emarginatus. They found higher variation in genetic diversity within populations than between populations. Here, we found that Sapindus was genetically diverse, with interspecific genetic variation significantly higher than intraspecific variation. S. mukorossi had higher levels of genetic variation and showed a pattern of geographic variation, whereas S. delavayi, S. rarak, and S. rarak var. velutinus showed low levels of intraspecific genetic variation. These results differ from those of Mahar et al., perhaps because their population originated in India, where the germplasm collection was smaller and unevenly distributed (Mahar et al., 2011a,b, 2013). It is also possible that interspecific incompatibility is lower between the Indian species. Overall, we believe that these results provide molecular biological support for the current consensus taxonomy of Sapindus in China. Furthermore, it suggests that special attention should be paid to Sapindus rarak var. velutinus with respect to its relationship with the other Sapindus species when carrying phylogenetic studies in future.

The construction of a core germplasm resource is an effective way to achieve efficient, scientific and rational conservation and utilization of genetic diversity (Xu et al., 2020), and the selection of an appropriate core population size is a crucial factor in establishing a core germplasm collection. Balakrishnan et al. (2000) and Zhang et al. (2010) suggested that the proportion of core germplasm should be determined by the size of the initial germplasm resource. Li et al. (2002) recommended sampling 5–40% of the core germplasm in crops, with 10% being optimal. However, woody plants are more genetically diverse, and the sampling percentage of core germplasm for woody plants is typically in the range of 10–45%. For instance, 14.71% (64/435) in apple (Malus domestica Borkh.) (Zhang et al., 2010), 17.96% (30/167) in Citrus reticulata (Garcia-Lor et al., 2017), 35% (247/707) in Eucalyptus cloeziana F. Muell (Lv et al., 2020), 35% (63/180) in Ginkgo biloba (Xuan et al., 2016), and 42.9% (300/700) in Chinese fir (Duan et al., 2017).

In this study, we selected a 40% (64/161) subgroup of the Sapindus core germplasm using Core Hunter. The population contained 53 S. mukorossi, nine S. delavayi, one S. rarak, and one S. rarak var. velutinus individuals (Supplementary Table 2). During the construction of core germplasm collections, allele retention is frequently considered as an evaluation criterion. For example, in Saccharum officinarum germplasm collection the criteria were to retain at least 70% of allele richness as well as other genetic diversity parameters (Balakrishnan et al., 2000). Compared with all Sapindus germplasm, the core germplasm collection had higher genetic diversity and maintained 80.66% of the allelic richness (Table 5), showing a balanced geographic composition (Supplementary Table 2). Previously, the core germplasm collection of lychee using 18 SSR markers resulted in 29.92% (38/127) individuals (Wang et al., 2012). The combined results indicate that the Sapindus core germplasm collection constructed in this study well represent the initial collection. The non-core germplasm of Sapindus population, also called reserve collection of Sapindus, is important for the conservation and utilization of Sapindus diversity, and it may harbor specific phenotypic, phenological and ecological characteristics to be of future use, thus it is important to conserve, exploit and understand the reserve collection as well. Although we have achieved our objectives of exploring genetic diversity and population structure in Sapindus and constructing a core germplasm collection using EST-SSR markers, our study still have some limitations. Compared with the millions of single nucleotide polymorphisms (SNP) and indel markers that can be obtained by whole genome resequencing or transcriptome sequencing, the number of EST-SSR markers in this study is relatively low, and the results obtained are not sufficient to fully elucidate the genetic structure of Sapindus. Hence, in future studies we hope to perform transcriptome sequencing or resequencing of this core germplasm collection to further explore interspecific and intraspecific genetic structure, evolutionary history, and regions under natural selection in Sapindus. Nevertheless, we believe the results of this study can help facilitate the efficient management, conservation, and utilization of Sapindus germplasm resources in the future.