To classify collected samples into distinct morphological groups, we initially used a subjective approach in accordance with a previous taxonomic study (i.e., Lu et al. 2020). According to studies conducted by Wang (2004) and Chao (2016), the population in east Taiwan that was not considered in the study by Lu et al. (2020) was regarded as A. daibuensis. We defined distinct morphological groups as a priori groups and used them in this study.

We determined whether a priori groups correspond to distinct taxa based on the concept of iterative taxonomy (Yeates et al. 2011), and species hypotheses were iteratively assessed using classical morphometrics and microsatellite genetic data. First, we performed morphometric analyses, including multivariate and univariate analyses, to evaluate quantitative morphological characters. Subsequently, we inferred genetic groups from microsatellite markers by using Structure software and explored genetic diversity and differentiation between genetic groups. Finally, we compared the results of the two datasets and concluded species delimitation.

For morphological comparison, living materials of Aspidistra were collected from Taiwan (Table 1, Fig. 1). A portion of these samples was converted into herbarium vouchers and deposited at TNU (The herbarium of National Taiwan Normal University, Taipei, Taiwan). For genotyping, the leaf material of these samples was dried using silica gel to preserve DNA and was stored at the Department of Life Sciences, National Taiwan Normal University, Taipei, Taiwan (Table 1, Fig. 1). In addition, herbarium specimens obtained from many herbaria (HAST, PPI, TAI, TAIF, and TNU; the acronyms are based on Index Herbariorum (https://sweetgum.nybg.org/science/ih/)) were examined directly, or digital plant specimens were evaluated online.

Figure 1. 

Sample collection location of the five a priori Aspidistra groups on the basis of morphological and geographical criteria. Color symbols: morphological and genetic data; white spots: only genetic data. Red triangle: AA, yellow inverted triangle: ADS, green square: ADE, fill diamond: AL, blue circle: AM.

In the first step, we selected 128 fresh specimens representing different a priori groups, including available flowering specimens (N = 97; Table 1, Fig. 1). We used 5 leaf characters and 16 floral characters for 97 samples (Tables 2, 3, Fig. 2). Raw data were normalized before analysis. Subsequently, we used a clustering method (clustering analysis [CA]) and an ordination method (discriminant analysis [DA]) to project and visualize trends for morphological variability across our samples, including leaf and floral characters. Finally, we determined whether the mean of quantitative traits significantly differed between morphological groups by using one-way analysis of variance (ANOVA). For characters determined to be significant, we again tested each pair of morphological groups through Tukey’s pairwise test and used letters to identify groups that differed significantly. All statistical analyses were performed using PAST statistical software v.4 (Hammer et al. 2001).

Figure 2. 

Morphological characters measured in this study A leaf characters B, C flower characters.

DNA was extracted from 15 to 25 mg of leaf material (stored in silica gel) for 128 samples by using the Viogene Plant Genomic DNA Extraction Miniprep System (Viogene-Biotek, Taiwan). On the basis of the study conducted by Lu et al. (2020), five populations without morphological data, namely H, HWS, NHL, DJL, and DLJ, were temporally classified into different a priori groups (Table 1). Nine microsatellite markers were amplified using the method of Huang et al. (2014). Genotyping was performed in a 48-capillary sequencer (Labnet MultiGene Gradient, Labnet International Inc., USA) by using 1.5 μL of DNA, 1000 μL of HiDi formamide (Thermo Fisher Scientific, Taiwan), and 8.3 μL of GeneScan-500 LIZ Size Standard (Applied Biosystems, Warrington, U.K.). To identify alleles (sizes of amplified PCR products), the resulting chromatograms were interpreted using Peak Scanner software v.1.0 (Applied Biosystems, Foster City, CA).

Genetic clusters were identified using the Bayesian clustering algorithm implemented in Structure v.2.3.4 (Pritchard et al. 2000) without a priori grouping. Ten independent runs were performed for each value of K ranging from 1 to 10 under a model assuming admixture and correlated allele frequencies (Falush et al. 2007). Each run comprised a burn-in period of 100,000 replications, followed by a run length of 1,000,000 Markov Chain Monte Carlo iterations. The results of replicated runs for each value of K from 1 to 10 were combined using Structure Harvester v.0.9.94 (Earl and von Holdt 2012). The optimal value of K was determined by calculating log-likelihood values and by using the ΔK method developed by Evanno et al. (2005). To visualize cluster assignments, the outputs of replicated runs were combined using CLUMPP v.1.1.2 (Jakobsson and Rosenberg 2007).

We used Arlequin v3.5.2 (Excoffier and Lischer 2010) to analyze molecular variance (AMOVA) and estimate the genetic diversity of each genetic group on the basis of the number of polymorphic loci as well as observed and expected heterozygosities (H_O and H_E) and the inbreeding coefficient (F). Differences in genetic diversity between genetic groups were characterized by computing the pairwise F_ST.

We identified five a priori groups on the basis of their morphological and geographical data. The names of these a priori groups were derived from one of the four recognized taxon names after matching with one (Table 1). (1) A total of 46 samples obtained from the higher mountain altitude of central to southern Taiwan were called “AA” because they corresponded to the species A. attenuata. (2) A total of 22 samples from central Taiwan (Defulan Trail and Shaolai Trail, Taichung City) were called “AM” because they corresponded to the species A. mushaensis. (3) A total of 14 samples from the Hengchun Peninsula were called “ADS” because they corresponded to the species A. daibuensis sensu Lu et al. (2020). (4) A total of 13 samples, which had distinctly thick perianth lobes, corresponded to the taxon A. daibuensis sensu Wang (2004) and Chao (2016) and were called “ADE.” (5) Finally, 2 samples from Huisun Experimental Forest Station, Nantou County, which had elongated stamen connectives, corresponded to the taxon A. longiconnectiva and were called “AL.”

Fig. 3A presents the results of the CA based on the Ward’s method that was performed to evaluate the 5 leaf characters and 16 floral characters for the 97 specimens. The morphological clustering groups almost fitted a priori groups, except for some mismatched samples. In the CA dendrogram, the samples were first divided into two groups: the first group consisted of samples from the AA group, and the second group consisted of samples from the ADE, ADS, AL, and AM groups. The second group was divided into four subgroups, which corresponded to the four a priori groups. However, except for the AL group, the other three groups included some mismatched samples, such as two ADE samples, and two ADS samples were assigned to the AM group (Fig. 3A). Axes 1, 2, and 3 of DA together accounted for 94.96% of the total variance. Axis 1 (59.81% of the relative contribution) was mainly determined on the basis of the ratio of the pistil height to the stamen height (PH/SH) and the distance from the perianth tube lobe base to the stigma apex (LBST). Axis 2 (20.9% of the relative contribution) was mainly determined on the basis of the distance from the perianth tube lobe base to the stigma apex (LBST), the ratio of the perianth tube length to the width of the perianth tube base (PTL2/PTBW), and the curve length of the stigma surface (SCL). Axis 3 (14.25% of relative contribution) was determined on the basis of the ratio of the blade length to the blade width (BL/BW), PH/SH, and the ratio of the leaf length to the blade length (LL/BL). In DA scatterplots, the AA and AL groups were unambiguously separated along axis 1, whereas the other groups were less pronounced (Fig. 3B, C). The ADE group was separated from the AM and ADS groups along axis 2, but the AM and ADS groups overlapped in the axis 1–axis 2 scatterplot (Fig. 3B). In the axis 1–axis 3 scatterplot, the ADS group could be distinguished from the AM group along axis 3, whereas the ADE group was located between these two groups (Fig. 3C).

Figure 3. 

Results of morphometric analyses A dendrogram of clustering analysis B, C scatterplot of discriminant analysis B scatterplot by Axis 1× Axis 2 C scatterplot by Axis 1× Axis 3.

A comparison of paired means by using one-way ANOVA revealed that the quantitative variables significantly differed among the five a priori groups, except for PTW2/PTBW and PTW3/PTBW (Table 3). The results of Tukey’s pairwise test for each pair of morphological groups were as follows. The AA group significantly differed from the ADE, ADS, AL, and AM groups in terms of three of the five leaf characters (LL, BL/BWL, and BL/BW) and the stigma disk angle (SDA). The AL group significantly differed from all the other groups in terms of floral characters (PTL2/PTBW and ratio of the perianth tube length [excluding the lobe] to the pistil height [PTL1/PH], and PTL1/SH). The leaf and floral characters differed among the ADE, ADS, and AM groups. In terms of leaf characters, leaf length (LL) significantly differed between the ADE group and the other two groups, BL/BW significantly differed between the ADS group and the other two groups, and the ratio of the leaf length to the petiole length (LL/PL) significantly differed between the AM and the other two groups. In terms of floral characters, the ADE had significantly thicker perianth lobes. In addition, the SCL, SDA, and stigma width significantly differed between the ADE group and the other two groups. Compared with the other two groups, the ADS group had a significantly smaller lobe base width, a smaller PTL1/PH ratio, and a PH/SH ratio.

The results of Bayesian structuring analysis performed using nine nuclear microsatellites in Structure software (Pritchard et al. 2000) revealed that the likelihood of data substantially increased with the number of imposed clusters K until K reached a value of 2 and a plateau was reached for a larger K value (Fig. 4A). Application of the delta K method developed by Evanno et al. (2005) demonstrated that the K value of 2 was the most likely number of clusters (Fig. 4B). We found that genetic cluster 1 matched the AA and ADS groups, and that genetic cluster 2 matched the AM, ADE, and AL groups. However, when investigating clustering solutions obtained at higher K values, a good match between the genetic clusters and a priori groups was observed at the K value of 4 (Fig. 4C). Cluster 1 matched the AA group, cluster 2 matched the AM and AL groups, cluster 3 matched the ADS group, and cluster 4 matched the ADE group. Therefore, we used the K value of 4 to analyze differences in genetic diversity between different genetic groups.

Figure 4. 

Results of genetic grouping A, B variation in the likelihood of data, L(K) (A) and of Delta K (B) as a function of the number of genetic groups (K) identified in 128 Aspidistra samples by using the software STRUCTURE C histogram of the genetic assignment of 128 Aspidistra individuals at K = 2 and K = 4. Individuals were ordered along the horizontal axis in accordance with their в group. Vertical bars indicate how the genome of each individual is partitioned into four clusters.

Furthermore, the genetic clusters indicated that the five populations without morphological data, namely H, HWS, NHL, DJL and DLJ, were temporally classified into different a priori groups according to Lu et al. (2020) have slightly different results. The DJL population was formerly classified to the ADS group, but now was allocated to the ADE group.

The genetic diversity of the genetic groups was substantial: the AL and AM groups displayed the highest diversity (H_E = 0.4403), whereas the ADE group exhibited the lowest diversity (H_E = 0.2842). The AA and ADS groups displayed intermediate diversity (H_E ranged from 0.3428 to 0.3951; Table 4A). The H_O of each genetic group was smaller than its H_E.

The inbreeding coefficient is a key parameter for understanding the number of matings between related individuals in a population. If no heterozygotes are present in a population, the inbreeding coefficient is 1.0. When the heterozygote frequency is equal to the Hardy–Weinberg expected value, the inbreeding coefficient is 0. Our results revealed that the inbreeding coefficients of the four genetic groups were significantly larger than 0, ranging from 0.1442 to 0.3705. This finding indicated that the degree of inbreeding varied within each group.

Pairwise F_ST is considered a satisfactory index to explain genetic differences among populations. Index values of 0–0.05, 0.05–0.15, 0.15–0.25, and >0.25 indicate nearly no, moderate, high, and significant genetic differences between populations, respectively (Wright 1978; Hartl and Clark 1997; Balloux and Lugon-Moulin 2002). Our results revealed that the pairwise F_ST between the genetic groups ranged from 0.2450 to 0.4939 (Table 4B). The pairwise F_ST of five of the six group pairs was >0.25. This result indicated the presence of significant genetic differences among the groups.

We determined differences in genetic diversity among the five a priori groups. The results were similar to the aforementioned results (Table 5). Except for the AL group, the H_O of the other four groups were all smaller than the H_E, and the inbreeding coefficient values were all greater than 0, ranging from 0.1442 to 0.3705. The pairwise F_ST of almost all the pairs was >0.25.

The CA and DA of leaf and floral characters revealed that the five a priori groups were well-supported by morphological groups (Fig. 3). Both the CA dendrograms and DA scatterplots demonstrated that the AA group differed from the AL group. The AA group corresponded to A. attenuata, and the AL group corresponded to A. longiconnectiva. With its longer leaves and concave stigma, A. attenuata exhibits clear morphological differentiation. Several authors have recognized all these characters to differentiate species (Hayata 1912; Liu and Ying 1978; Yang et al. 2001; Wang 2004; Chao 2016; Lu et al. 2020). A. longiconnectiva is characterized by deep perianth lobes, relatively broad and shallow perianth tubes, and elongated stamen connectives (Lu et al. 2020). The ADE, ADS, and AM groups corresponded to the eastern population of A. daibuensis, the southern and southeastern populations of A. daibuensis, and A. mushaensis, respectively. However, some mismatched samples were noted within these three groups. These mismatched samples may be attributed to the ambiguous morphological characters between the two species. This result is consistent with those of previous taxonomic studies. Wang (2004) and Chao (2016) have considered that the southern and southeastern populations of A. daibuensis correspond to A. mushaensis, and that the eastern population of A. daibuensis corresponds to A. daibuensis. However, Lu et al. (2020) believed that A. mushaensis exists only in central Taiwan, and that the southern and southeastern populations of A. daibuensis correspond to A. daibuensis and not A. mushaensis. Whether the ADS group corresponds to A. daibuensis or A. mushaensis remains unclear. However, our morphometric analyses revealed that the ADS group differed from the AM and ADE groups in terms of a wider blade width, a smaller lobe base width, more stigma present near the perianth opening, and the stamen attached on the perianth tube base or near the base. Thus, we suggest that the ADS group should not be classified as the AM or ADE group.

The results of Bayesian structural analysis revealed that a K value of 2 was the most favorable for genetic clustering. Genetic cluster 1 consisted of the AA and ADS group, and genetic cluster 2 consisted of the ADE, AL, and AM groups. However, this genetic division was inconsistent with the morphological division based on CA, wherein in the CA dendrogram, the samples were divided into the AA group and the other groups (ADE, ADS, AL, and AM). Although two clusters were obtained in our dataset based on ∆K, we determined that a K value of 4 led to more favorable matching with the a priori groups. Demographic, environmental, and historical processes are multifaceted and complex and have led to different organization levels in the genetic structure of species. Therefore, Meirmans (2015) suggested that different K values may reflect different demographic processes, and that a biologically interpretable pattern obtained from a suboptimal K value is better than a completely unrealistic pattern obtained from the optimal K value. Thus, we used the K value of 4 to determine the relationship between genetic and morphological groups. Although genetic clusters almost supported our morphological groups, AL clustered with AM to form a genetic group. This result can be attributable to the small sample size and incomplete lineage sorting. First, the AL group corresponded to A. longiconnectiva, which is a rare species in Taiwan, and we only sampled six individuals from a locality. The Structure software usually fails to identify genetic groups represented by few individuals, even when one group is well-represented (Porras-Hurtado et al. 2013; Wang 2017). Second, AL and AM are morphologically similar (Lu et al. 2020) and are sympatrically distributed. We believe that AL and AM are closely related and may have diverged recently. Genetic data also indicated that they were closely related. Thus, the lineage sorting of these two species is incomplete.

The inbreeding coefficient indicates the effect of inbreeding on homozygosity by quantifying the deviation in observed genotypic frequencies from those expected under the Hardy–Weinberg equilibrium (Charlesworth 2003). Inbreeding, null alleles, and undetected genetic structure increase observed homozygosity and the inbreeding coefficient, whereas outbreeding, mutation, and inbreeding depression tend to reduce the inbreeding coefficient (Stoeckel et al. 2006; Waples 2018). Our results revealed that the inbreeding coefficients of all the genetic clusters were significantly greater than 0. Therefore, each genetic group had a high possibility of inbreeding. Inbreeding promotes population fragmentation and provides a selective advantage at the population level, resulting in structured populations, which contribute to speciation (Joly 2010; Olsen et al. 2020).

In our genetic differentiation analysis, we observed significantly high pairwise F_ST values, ranging from 0.2450 to 0.4939. Among the pairwise F_ST values of the six genetic group pairs, the AA with ADE pair and the AA with AL+AM pair exhibited high genetic differentiation (0.4939 and 0.4105, respectively), whereas the pairwise F_ST values of the AA with ADS pair was slightly smaller (0.2450) but still indicated high genetic differentiation (Wright 1978; Hartl and Clark 1997; Balloux and Lugon-Moulin 2002). The remaining genetic group pairs, namely ADE with ADS, ADE with AL+AM, and ADS with AL+AM, had significantly high pairwise F_ST values (0.3774, 0.2710, and 0.3761, respectively). We believe that these three groups have high genetic differentiation. This result also supports our morphometric analysis finding that ADE, ADS, and AM are distinct entities.

We performed an iterative analysis of morphometric and microsatellite markers to define, assess, and delineate species within the Aspidistra genus in Taiwan. The integration of the two analyses was almost congruent for the five a priori groups. In addition, the population genetic data indicated that frequent inbreeding and high genetic differentiation may shape the species diversity of Aspidistra in Taiwan. Therefore, we suggest that the five a priori groups should be regarded as five different taxa. Considering the morphological similarities between these five taxa, their geographical distribution as well as the preliminary ranking rule providing by Tillich and Averyanov (2018), we proposed the following taxonomic treatments to update the taxonomy of Aspidistra in Taiwan: AA corresponds to the species A. attenuata, which is distributed in the mountain range at higher altitudes in the central and southern Taiwan; AM corresponded to the species A. mushaensis, which is only distributed in central Taiwan; AL corresponded to the species A. longiconnectiva, because it is genetically similar to A. mushaensis, we reduced it to the variety level; ADE corresponded to the species A. daibuensis, which has restricted distribution in southeastern Taiwan; ADS morphologically resembles ADE, but they differ in their genetic data. Therefore, we considered regarding it as a new variety of A. daibuensis.