Monitoring the population-level genetic diversity for a species offers invaluable insights into conservation practices as the variation directly reflects the evolutionary and ecological process that the species has experienced (Olson et al., 2016; Mimura et al., 2017; Schierenbeck, 2017; de Kort et al., 2021). Biodiversity encompasses variations at all levels from the infraspecific to the community, yet infraspecific diversity, including population-level genetic variation, has been relatively underestimated despite its rate of loss exceeding that of higher-level diversity (Mimura et al., 2017; Kooyman et al., 2020). Within the last century, large bodies of land surfaces have been heavily exploited by human activities resulting in accelerated rates of biodiversity loss (Hautier et al., 2015; Scheffers et al., 2016). Accordingly, the role of population genetic diversity has become increasingly significant in terms of responding to biodiversity loss and subsequent conservation strategies (des Roches et al., 2021; Hohenlohe et al., 2021; Stange et al., 2021). For example, Lee et al. (2018) suggested risk in deleting a rare myrmecochorous herb, Plagiorhegma dubium from the Red List of South Korea based on the pattern of population genetic diversity of simple sequence repeat data despite the recent growth of census population sizes.

Orchidaceae is a cosmopolitan plant group that comprises well beyond 20,000 taxa exhibiting the most remarkable diversity of all levels from infraspecific to generic (Summerhayes and Withner, 1960; Cozzolino and Widmer, 2005). Ironically, at the same time, orchids are one of the most rapidly disappearing plant groups, which are listed heavily in the IUCN Red List, for example, more than 800 orchids are designated as threatened (Cribb et al., 2003; Wraith et al., 2020; IUCN, 2021). Coupled with intrinsic factors such as complex symbiotic relationships with specialized pollinators and mycorrhizal fungi, anthropogenic threats such as climate change, habitat destruction, and poaching expedite the rate of diversity loss in Orchidaceae (Swarts and Dixon, 2009; Phillips et al., 2020; Wraith et al., 2020). Given the current pace of climate changes, commonly required specific niche preferences of terrestrial orchids will likely expose them to high extinction vulnerability (Orchid Specialist Group [UCN/SSC], 1996; Swarts and Dixon, 2009). Local populations of terrestrial orchids are likely to become isolated and patchy under ongoing environmental changes. Consequently, the infraspecific genetic diversity would drop more quickly than the species and generic-levels diversities. Therefore, rigorous monitoring of population-level genetic diversity in terrestrial orchids is highly recommended for conservation and management plans.

Habenaria Willd. is one of the most diverse taxonomic groups of the terrestrial orchids distributed across the old and the new worlds, which consists of more than 800 taxa (Govaerts et al., 2010). Despite the large breadth of distribution, ∼10% of Habenaria species display a high vulnerability to extinction (IUCN, 2021). The proportion of threatened taxa in the genus will likely grow over time due to ongoing environmental changes and a lack of knowledge on many Habenaria taxa. Despite the increased vulnerability, there are only a handful of Habenaria taxa for which the population-level diversity pattern has been investigated (Zhang and Gao, 2017; Chung et al., 2018; Tachibana et al., 2021). Habenaria linearifolia (Maximowicz) Szlachetko, a rare terrestrial orchid with a narrow distribution (South Korea, E China, Japan, and Russian Far East) is restricted to bogs and wet grasslands (Lee and Choi, 2006; Chen and Cribb, 2009; Chung et al., 2018). As for many terrestrial orchids, populations of H. linearifolia are small, patchy, and isolated (Chung et al., 2018), which likely promotes extinction risk. Notably, the species is not on any protected species list, partly due to a lack of studies examining the infraspecific diversity status. Accordingly, an infraspecific genetic diversity study at a genomic scale might offer a diagnosis of the current level of diversity and vulnerability for H. linearifolia.

In this study, we surveyed genomic diversity and the spatial pattern of H. linearifolia, a rare terrestrial orchid in South Korea. Habenaria cruciformis, the most closely related species having a questionable taxonomic boundary problem with H. linearifolia was also included in the study. Using thousands of genome-wide single nucleotide polymorphism (SNP) markers, we aimed to (1) assess the extent of genetic diversity within and among populations; (2) evaluate gene flow among H. linearifolia populations and between the two species; (3) determine the presence of genetic bottlenecks; and finally (4) identify any cryptic diversity that is presently unknown. Given the specialized habitat preferences and the obligatory symbiotic relationships with pollinators and mycorrhizal fungi, we hypothesized that the populations of H. linearifolia have largely differentiated with low within-population genetic diversity. Habenaria cruciformis population might have largely diverged from all H. linearifolia populations with very limited gene flow if the species is an independent species with complete reproductive isolation. The results may suggest separate management units and/or assist in determining the conservation priority.

In total, the library produced ∼108 Gbp with an average GC content of 39%. From the initial SNP call, we isolated over 91K SNPs. After a series of pruning steps with the thresholds described, 2,734 SNPs were retained for the downstream analyses. Overall genetic diversity of H. linearifolia varied across populations although the differences were not significant (Table 1). The number of alleles adjusted by the population size ranged from 1.44 (GW and JJ; see Table 1 for the population acronyms) to 1.58 (AH and NW). The mean He (= 0.26) and percentage of polymorphic loci (mean = 70.9) were the lowest in HS at which the collected sample was the smallest (Table 1). Among-population level genetic divergence measured by F_ST varied largely across the population pairs ranging from 0.07 (JJ/SH) to 0.30 (HS/SH; Table 2). All F_ST values were statistically significant at the P < 0.01 level.

According to AMOVA, ∼60% of genetic variance was attributable to variations within individuals (Table 3). A fairly large proportion of total genetic variance (∼18%) was partitioned among regional groups redefined by the Bayesian assignment analysis (Table 3). Genetic variance among populations within each region was smaller than that derived from among regions (Table 3). We found a significant IBD pattern from the Mantel test result (r = 0.52; P < 0.01; Figure 1).

Our clustering results of both PCoA and STRUCTURE exhibited similar patterns, which largely separates three populations in Jeju Island from the remaining eight populations (Figures 2, 3). The PCoA plot of PC1 and PC2 divided the populations into four clusters primarily by the PC1 axis which explains about 13% of the total variance (Figure 2). Notably, PC1 separates the three Jeju populations of H. linearifolia into two clusters, whereas one H. cruciformis population represented as YW was isolated from the rest of the populations exclusively by the PC2 axis (∼8%; Figure 2). The best K determined by ΔK was K = 2 (Supplementary Figure 1); however, we plotted K = 2–9 and presented the bar charts here to explore the admixture patterns as the number of K clusters increased (Figure 3). At K = 2, three Jeju populations formed a regional deme sharing genotype frequencies represented mostly by an orange ancestry with a blue ancestry admixed (Figure 3A). As the number of the ancestry groups grew to three (K = 3), H. cruciformis started exhibiting a unique ancestry pattern with a minor portion of admixture (Figure 3B). When the K increased to nine, the pattern of the shared ancestry reflected the geography such that H. linearifolia populations formed four regional demes, that is, Jeonnam, Gyeongnam, Jeju, and Gangwon (Figures 3C,D).

The indirect measures of gene flow estimated from Migrate-N greatly varied among the demes and were somewhat consistent with the clustering patterns (Table 4). Overall, the gene flow among the regional demes showed asymmetry, in particular, it was the strongest between demes 1 (JN) and 5 (GW) (5 → 1 M = 0.258; 1 → 5 M = 0.014; Table 4). Deme 5 (GW; see Table 1 for the deme information) appeared to migrate much more to the other demes while it received a very limited number of immigrants from the rest of the demes (Table 4). However, the statistical support of the rates of emigration to the other demes from Deme 5 was weak (Table 4). Deme 1, Deme 2, and Deme 3 showed a rather high migration rate into the remaining demes (∼0.2–0.3/generation/site; Table 4). Based on the Bottleneck results, all populations have experienced significant reductions in population sizes with-in the recent past (Table 5).

Terrestrial orchids, one of the species-rich groups, are rapidly disappearing from the earth due to their specialized habitat preference and their complex symbiotic relationships with mycorrhizal fungi and pollinators (Swarts and Dixon, 2009; Phillips et al., 2020; Wraith et al., 2020). Under the rapid and ongoing biodiversity loss induced by anthropogenic environmental changes, the reduction in orchid diversity is likely accelerated. By monitoring the population-level genetic diversity of two rare terrestrial orchids, we revealed three major findings. First, the orchid populations have experienced genetic bottlenecks and were genetically isolated with a rather low gene flow, which might suggest the possibility of local extinction (Tables 2, 4). Despite the limited gene flow and reductions in population sizes, the within-population genetic diversity in H. linearifolia is comparable to the other terrestrial orchids (Table 1). Finally, we found cryptic variation in the regional populations on Jeju Island from the Bayesian assignment analysis that was not shared with the inland populations (Figure 3).

Overall, the diversity of 2,734 SNPs observed in H. linearifolia and H. cruciformis appeared to be comparable or even greater to that estimated for other terrestrial orchids [Table 1; Pleione formosana, Ho = 0.27 (Chao et al., 2021); and Epipactis helleborine Na = 1.21, He = 0.054 (Jacquemyn et al., 2020)]. The amount of genetic variability harbored in a rare orchid such as H. linearifolia was surprising considering the population bottlenecks we found (Table 5). In addition, it is noticeable that the high molecular variation was observed only in the degree of heterozygosity and not in allelic richness (Table 1). Populations that have experienced a sudden demographic crash are expected to be low in genetic variation due to the massive loss of alleles resulting from genetic drift, yet the decline of genetic diversity may not be prominent if the bottleneck was a single event (Nei et al., 1975; Carson, 1990). Given the lengthy longevity of orchids, the one migrant per generation rule (Mills and Allendorf, 1996) might alternatively be the causal mechanism for the maintenance of diversity. However, migration rates lower than one (mean migration rate = 0.14; Table 4) undermines the probability of the hypothesis. The high heterozygosity observed might also be resulted from hybridization among multiple lineages. Indeed, Thompson et al. (2010) stressed the importance of hybridization as a source of genetic variation for rare plant species. However, our data were not equipped to test the possibility as we could not sample all regional lineages due to the inaccessibility of some northern and island populations. Further study may expand the sample sizes and determine the true mechanism for the maintenance of genetic variability.

Pairwise genetic divergences values (mean pairwise F_ST ≈ 0.2) were greater than the ones observed in many rare orchids [Spiranthes spiralis, mean F_ST = 0.02 (Machon et al., 2003); Caladenia huegelii, F_ST = 0.047 (Swarts et al., 2009); Orchis purpurea, F_ST = 0.09 (Jacquemyn et al., 2020)], which clearly indicate large population divergence of H. linearifolia (Table 2). Although reduced gene flow has been speculated based on the rarity and narrow distribution of most orchids, a recent meta-analysis revealed rather contrasting patterns of gene flow (Phillips et al., 2012). Considering the capability of long-distance dispersal resulting from the minute size and massive quantity of seeds in orchids, moderate-to-high gene flow found in the meta-analysis indeed is not surprising. Habenaria orchids also produce thousands of dust-like seeds (Chen and Cribb, 2009), thus the large population differentiation observed here for H. linearifolia was unexpected. The increased population divergence might be attributed to the limited gene flow we observed in South Korean Habenaria (Table 4). The minute orchid seeds may travel far, but the actual movement distance will vary for each seed produced and is highly dependent on the neighboring vegetation height (Jersáková and Malinová, 2007; Brzosko et al., 2017). As South Korean Habenaria mostly resides in grasslands surrounded by tall grasses, low gene flow limited by seed dispersal seems likely. Alternatively, such population diversification could also arise by a combination of strong genetic drift and natural selection (Gentry and Dodson, 1987; Tremblay et al., 2005). The effect of genetic drift is expected to some extent, however, with the current data set, it is uncertain how natural selection could impact the process of population diversification. To tease apart the role of limited gene flow from the contribution of drift and selection on population diversification, empirical studies assessing selection such as common garden experiments are needed.

Notably, Bayesian clustering analysis revealed an unexpected assignment pattern. Based on the best K (K = 2), 2,734 SNPs have the eleven populations of H. linearifolia and H. cruciformis assigned to two demes consisting of the inland deme and the Jeju Island deme (Figure 3A). It is rather surprising that H. cruciformis represented by the Youngwol population was not separated from the inland populations of H. linearifolia. As the cluster number, K, grows by 3, H. cruciformis started being divided from H. linearifolia showing a distinctive pattern of genotype frequency. However, it is noticeable that Jeju populations also exhibited a rather unique pattern of genotype frequency with some level of admixture in K3 (Figure 3B). Habenaria cruciformis was first reported by Ohwi (1951) based on the morphological aforementioned characters (the spur length and the shape of lateral lobes). On the other hand, Jeju populations were never recognized as independent taxon due to a lack of measurable morphological differences. Except for the Jeju populations, most H. linearifolia populations diversified following a significant IBD pattern (Figure 1).

Due to contemporary and past evolutionary processes such as long-term population isolation, single species might diverge into multiple evolutionarily significant units [ESUs; refer (Fraser and Bernatchez, 2001) for the terminology details] or even harbor a cryptic species (Bickford et al., 2007; Warner et al., 2015; Hill et al., 2016; Martin et al., 2016). With the recent advancements of modern molecular techniques, signatures of ESUs and cryptic species have been mounting. Despite the high morphological similarities, our molecular data exhibited clear separation of the island deme from both inland H. linearifolia and H. cruciformis suggesting a possibility of cryptic ESU or even a species. Delimitation of the formerly undescribed ESUs and/or cryptic species provides critical insights by identifying genetically cohesive units. To confirm the emergence of a cryptic taxon, a rigorous examination of morphological characters and a phylogenetic verification should be added. Although cautions should be taken, our data highly suggest that the Jeju Island deme must be considered as an important ESU.

The genetic diversity pattern we observed indicated that the species has been experiencing some level of diversity loss and isolation among populations with low gene flow. Habenaria linearifolia requires a specific type of habitat (wet meadows and bogs) that is currently shrinking with the anthropogenic interferences (Ohwi, 1965; Essl et al., 2012). Without appropriate conservation strategies, the species distribution and size of each population are likely going to diminish even further. For the conservation practices, we suggest building both close monitoring and protection programs, particularly for the newly discovered ESUs. Preparation for ex situ conservation, such as collecting and preserving seeds, is also recommended along with intense ecological monitoring.

The original contributions presented in the study are publicly available. This data can be found here: National Center for Biotechnology Information (NCBI) BioProject database under accession number PRJNA758713.

S-RL and S-YJ conceived the research idea and designed the experiments. S-YJ got the funding and revised the manuscript. T-YC conducted library preparation and related laboratory work. S-RL and T-YC performed genetic and statistical analyses. S-RL wrote the manuscript. All authors contributed to the article and approved the submitted version.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.