﻿Aspidistradaibuensis var. longkiauensis, a new variety of Aspidistra (Asparagaceae) from Taiwan, identified through morphological and genetic analyses

﻿Abstract Aspidistra Ker Gawl. is one of the the most diverse and fastest-growing genera of angiosperm. Most Aspidistra species have been discovered in a limited area or a single site through morphological comparison. Because of the lack of population studies, morphological variation within species and the boundaries of some species remain unclear. In recent years, combining genetic and morphological markers has become a powerful approach for species delimitation. In this study, we performed population sampling and integrated morphometrics and microsatellite genetic diversity analyses to determine the species diversity of Aspidistra in Taiwan. We identified three species, namely Aspidistraattenuata Hayata; A.daibuensisHayatavar.daibuensis; A.mushaensisHayatavar.mushaensis; and reduced A.longiconnectiva C.T.Lu, K.C.Chuang & J.C.Wang to the variety level, and described a new variety, A.daibuensisHayatavar.longkiauensis. The description, diagnosis, distribution, and photographs of this new variety as well as a key to the known Taiwanese Aspidistra are provided.


Introduction
The genus Aspidistra Ker Gawl. (Asparagaceae) is native to eastern and southeastern Asia, particularly China and Vietnam (Tillich 2014), and Aspidistra species typically grow under forest canopies and shrubs in high-rainfall areas. Plants of Aspidistra are characterized by a perennial herbaceous habit, conspicuous rhizomes, solitary or two to four tufts' leaves, a variety of fruits and a highly diversified flower structure (Li et al. 2000;Lin et al. 2010). Aspidistra is one of the fastest growing genera in Angiosperm (Tillich and Averyanov 2018). In the past 15 years, the number of Aspidistra has considerably increased from approximately 93 taxa (Tillich 2008) to more than 200 species (Kalyuzhny et al. 2022). So far, new species of this genus have been discovered continuously. However, phylogenetic analysis of this genus has not yet been complete performed.
In Taiwan, Hayata (1912Hayata ( , 1920 described three species, namely Aspidistra attenuata Hayata, A. daibuensis Hayata, and A. mushaensis Hayata. However, the boundaries of some of these species are inconsistent. For example, Ying (2000) described that these three species should be combined into A. elatior Blume and regarded them as a variety, namely A. elatior var. attenuata (Hayata) Ying. Subsequently, Yang et al. (2001) demonstrated that the previously described three species differed from A. elatior, and that A. mushaensis is a synonym of A. attenuata Hayata. Wang (2004) and Chao (2016) have accepted the species concept of Hayata and reported that A. daibuensis is distributed in south and east Taiwan. However, Lu et al. (2020) reported that A. daibuensis is restricted to only south Taiwan and did not consider the population in east Taiwan. Wang (2004) suggested that A. mushaensis is distributed in not only central Taiwan but also south Taiwan. However, Lu et al. (2020) disagreed with this suggestion. A recent study suggested that the genus Aspidistra contains four species: A. attenuata Hayata; A. daibuensis Hayata; A. longiconnectiva C.T.Lu, K.C.Chuang, & J.C.Wang; and A. mushaensis Hayata (Lu et al. 2020).
Taxonomists have recognized nearly all Aspidistra species through the morphological comparison. However, this approach of species delimitation can cause some identification problems because of the intraspecific morphological variation or small morphological differences between closely related species (Pratt and Clark 2001;Whittall et al. 2004). Advances in many fields, such as molecular genetics, have helped taxonomists determine species boundaries and identify cryptic species (Whittall et al. 2004;Ellis et al. 2006;Bickford et al. 2007).
Population genetic methods, such as the algorithm developed by Pritchard et al. (2000) that was implemented in Structure software, can be used to identify welldifferentiated genetic clusters and thus detect gene flow barriers. These methods can be employed to detect distinct species even if they are not yet reciprocally monophyletic at many genes due to incomplete lineage sorting (Duminil et al. 2015). In addition, by determining associations based on morphological data and considering the spatial distribution of detected genetic groups, we can investigate whether these genetic groups correspond to distinct species. Recently, many studies, such as those on Ancistrocladus Wall. (Turini et al. 2014), Asteropyrum J. R. Drumm. et Hutch. (Cheng et al. 2021), Greenwayodendron Verdc. (Lissambou et al. 2019), Polygonum L. (Mosaferi et al. 2016), and Santiria Blume (Ikabanga et al. 2017), have reported that the combination of morphological and population genetic data can be used to delineate species complexes.
In this study, we integrated morphological and genetic data to evaluate species delimitation within the genus Aspidistra in Taiwan. In addition, we conducted population genetic and morphometric analyses to confirm the previous classification hypotheses and to determine the number of species in the genus Aspidistra in Taiwan. We further performed genetic differentiation analysis for groups to determine whether they are the conspecific populations or distinct species.

Recognition of a priori species groups
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.

Plant material sampling
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.

Identification of morphological groups
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).  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).  The ratio of the perianth tube length to the width of the perianth tube base PTW1/PTBW The ratio of the perianth tube width on the stamen-attached position to the width of the perianth tube base PTW2/PTBW The ratio of the perianth tube width on the stigma apex to the width of the perianth tube base PTW3/PTBW The ratio of the perianth tube width on the lobe base to the width of the perianth tube base PTL1/PH The ratio of the perianth tube length (exclude lobe) to pistil height PTL1/SH The ratio of the perianth tube length (exclude lobe) to the height of the stamen attached on the perianth tube (stamen height) PTL1/LBL The ratio of the perianth tube length (exclude lobe) to lobe length PTL1/PTL2 The ratio of the perianth tube length (exclude lobe) to the perianth tube length PH/SH The ratio of pistil height to the height of the stamen attached on the perianth tube (stamen height) The letters after the numbers were used to identify groups with significantly differences.

Identification of genetic groups
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).

Genetic diversity and genetic differentiation between genetic groups
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 .

Definition of a priori groups
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." Morphological differentiation between a priori groups 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). 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.

Identification of genetic groups
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.
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.

Genetic diversity and genetic differentiation between genetic groups
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.

Comparison with morphological and a priori groups
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.

Comparison with genetic groups and a priori groups
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.

Genetic diversity and genetic differentiation between genetic groups
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.

Integration of morphometric and genetic data
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.
Etymology. The specific epithet "longkiau" means that the geographical distribution of this species is mainly distributed in Hengchun Peninsula, Pingtung County, Taiwan. The area was known as "Longkiau" in early records.
Phenology. Flowering from February to June. Distribution and habitat. Aspidistra daibuensis var. longkiauensis is geographically confined to the Hengchun Peninsula, the most southern part of Taiwan and the southeastern part of Taiwan. The population typically grows on slopes under the canopies of primary or secondary forests. It is distributed at an elevation of 200-500 m.
Conservation status. Aspidistra daibuensis var. longkiauensis is known from three populations in Hengchun Peninsula, Pingtung County, Taiwan. Based on the specimen records, the area of occupancy (AOO) is ca. 15 km 2 by GeoCAT (Bachman et al. 2011). Following the criteria of IUCN (2019), this species is assessed as endangered (EN B2ab(ii)).
Notes. This new variety is regarded as A. daibuensis (Chang and Hsu 1974) or A. mushaensis (Wang 2004;Chao 2016). We determined that this species does not have thick perianth lobes as A. daibuensis var. daibuensis (Hayata 1912). In addition, A. daibuensis var. longkiauensis differs from A. mushaensis var. mushaensis in terms of the ratio of the blade length to the blade width, perianth lobe base width, the distance from the lobe base to the stigma apex, and the ratio of the pistil length to the stamen height ( Table 6). The stamen of this new species is inserted on the perianth tube base or near the base instead of the low part of the perianth tube, as observed in the previous former two species. We considered it can be distinguished from A. daibuensis var. daibuensis and A. mushaensis var. mushaensis. However, it is morphologically more similar to A. daibuensis, so we considered regard it as a variety of A. daibuensis.
Additional specimen examined. Taiwan  Leaf length to leaf blade length < 1.2; the perianth tube length to the width of the perianth tube base < 0.9; the perianth tube length (exclude lobe) to the height of the stamen attached on the perianth tube < 2.