The target species is the root-hemiparasitic herb species, P. hallaisanensis ( Figures 1C, D ). This species is endemic to Korea and is legally protected as an endangered species (grade II) of the Ministry of Environment (Chung et al., 2023). It originally inhabited Korean mountaintop grasslands (altitude: 1,400–1,500 m) such as Hallasan and Gayasan (Kim et al., 2018), with voucher specimens stored at the herbariums of Inha University (Cho. 98454) and the Korea National Park Research Institute (Gaya_20160439), respectively (Cho and Choi, 2011; Han et al., 2022). However, the most wild populations have disappeared due to habitat loss and climate change (Kim et al., 2018). This species is known to be phylogenetically close to other hemiparasitic Orobanchaceae species, such as Pedicularis spicata, Pedicularis verticillata, Pedicularis alaschanica, and Pedicularis longiflora, according to previous ribosomal and chloroplast DNA studies (Cho and Choi, 2011; Cho et al., 2018; Zhang et al., 2020). Morphologically, P. hallaisanensis is distinguishable by its dense glandular hairs covering the entire body, a shorter galea compared to the lower lip, and three unequal calyx lobes when compared to other allied Orobanchaceae species (Cho and Choi, 2011; Kim et al., 2018). Although P. hallaisanensis has traditionally been speculated to be either an annual, biennial, or perennial species (Cho et al., 2018; Kim et al., 2018), our monitoring was the first to confirm that this species features a strict biennial lifecycle (see our Results section).

The study area was located in a mountaintop grassland within Gayasan National Park in South Korea (35°49′25″N, 128°7′10″E) ( Figures 1A, B ). The altitude of the study area ranges from 1,410 to 1,430 m above sea level (asl), with slopes varying between 0° and 20°. The average annual precipitation is 1,296 mm, with a relative humidity of 74.8% and an average air temperature of 7.6°C. The soil is less than 20 cm, lacking distinct evidence of horizon development (entisols), and is underlain by an impermeable bedrock layer. This was the only study area available for studying the target species, although we have investigated all known natural habitats around Hallasan, Gayasan, Seoraksan, and Bangtaesan from 2019 to 2023 in search of additional wild P. hallaisanensis populations. National inventory data from the National Institute of Ecology have also recorded no wild P. hallaisanensis population remaining in any natural habitats since 2020, except for this study area.

All P. hallaisanensis individuals were counted in August from 2021 to 2023 to quantify the population size in the study area. Of these, one individual was used for the draft genome analysis, and 20 were used for SNP analyses. All P. hallaisanensis individuals in the study area were sorted into age classes according to the elapsed time after germination (first-year seedling and second-year adult). Each P. hallaisanensis was labeled and monthly monitored to track any morphological changes throughout the biennial lifecycle during the growing season (April–November).

Fresh leaves of a P. hallaisanensis individual were sampled in June 2023, stored in an icebox (4°C), and brought to the laboratory for draft genome assembly. Genomic DNA was extracted from the leaf samples using the Aprep Total DNA KIT (APBIO, Namyangju, South Korea) based on the manufacturer’s instructions. The extracted DNA was then quantified with a Thermo Scientific Nanodrop 8000 spectrophotometer (Fisher Scientific, Waltham, MA, USA), digested using the ApeKI enzyme (GCWGC), and reorganized into short reads of 151 bp length, following the protocols of Elshire et al. (2011) and Oh et al. (2023). Sequencing was conducted using the Illumina Hiseq X Ten Platform (Illumina Inc., San Diego, CA, USA). Illumina raw reads were filtered using Trimmomatic v.0.39 to exclude poor-quality reads (window size: 4, mean quality: ≧ 15, leading and trailing: ≧ 3, read length: ≧ 36 bp) (Bolger et al., 2014), and de novo assembly was done using SOAPdenovo2 v.2.04 (K-mer = 69) (Luo et al., 2012).

Two 5 × 5 m plots were established for SNP analyses in two different microhabitats (MH-1 and MH-2). Although the distance between these two microhabitats was only 20 m, they differed in terms of historical anthropogenic interventions ( Figure 1B ). MH-1 had been artificially flattened and managed as a heliport until the early 2010s ( Supplementary Figure S1 ), but it has recently become revegetated after such management was stopped. Conversely, MH-2 was located outside the old heliport sites and was relatively sheltered from heavy anthropogenic interventions. Both MH-1 and MH-2 contained both odd-year-flowering (OYF) and even-year-flowering (EYF) P. hallaisanensis cohorts, in contrast to several other microhabitats, which included either OYF or EYF cohorts only.

Fresh cauline leaves from five OYF and five EYF were randomly sampled from each plot for SNP detection in June 2023 and 2022, respectively (n = 20). Short DNA reads (151 bp) were obtained using the protocols of Elshire et al. (2011) and Oh et al. (2023) and sequenced with the Illumina Hiseq X Ten platform (Illumina Inc., CA, USA). Low-quality raw reads were then removed using cutadapt v.1.8.3 (Li, 2013) and Trimmomatic v.0.39 (Bolger et al., 2014). Filtered clean reads were mapped to the assembled draft genome of P. hallaisanensis using BWA v.0.7.17-r1188 (Li, 2013), and raw SNPs were detected using SAMtools v.0.1.16 (Li et al., 2009). To ensure SNP quality, only biallelic SNP loci without any missing values throughout all 20 samples were selected for further statistical analyses (3716 SNPs in total).

To describe the genetic diversity, Nei’s genetic diversity (GD), polymorphism informative content (PIC), minor allele frequency (MAF), and observed heterozygosity (Ho) were calculated using the snpReady package in R v.4.3.2 software (Granto et al., 2018).

Permutational multivariate analysis of variance (PERMANOVA) and permutational analysis of multivariate dispersion (PERMDISP) were conducted using Bray–Curtis dissimilarity based on 9,999 permutations to test the effects of flowering group (OYF or EYF) and microhabitat (MH-1 or MH-2) on the multivariate genetic centroids and dispersions of the SNP data from 20 P. hallaisanensis samples (α = 0.05). Nonmetric multidimensional ordination scaling (NMDS) was further conducted using Bray–Curtis dissimilarity to visualize the multivariate genetic variability, and a general linear model (GLM) was applied to test the relationship between NMDS axes, flowering group, and microhabitat (n = 20, α = 0.05). These analyses were performed using the vegan (Oksanen et al., 2024) and agricolae (De Mendiburu and Simon, 2015) packages in R v.4.3.2 software (R Core Team, 2023).

Neighbor-joining method and k-means clustering of the silhouette width approach were implemented using the factoextra (Kassambara and Mundt, 2016) and ape (Paradis et al., 2004) packages in R v.4.3.2 software to show genetic differences among the sampled P. hallaisanensis (n = 20) (Batool and Henning, 2021). STRUCTURE 2.3.4 software was also used with 10,000 burn-in periods, 100,000 Markov chain Monte Carlo (MCMC) replications, and 10 iterations to estimate the membership probability of each P. hallaisanensis sample according to hypothetical ancestral genotypes (Pritchard et al., 2000).

Our monitoring demonstrated that P. hallaisanensis strictly required 2 years of lifecycle per generation, with approximately 8 and 10 months of growing periods as first-year seedling and second-year adult, respectively. From April to May, first-year seedlings of P. hallaisanensis germinated and developed leaves and roots, eventually forming overwintering buds in November ( Figure 2A ). There were only rosette leaves without distinguishable shoots and cauline leaves in the aboveground of the first-year seedlings ( Figure 2A ). However, second-year adults of P. hallaisanensis rapidly established shoots and cauline leaves from the overwintering buds starting in April. Flowering and seed production of occurred from August to November in the second year, and all second-year adults died after producing seeds ( Figure 2A ). Facultative annual and perennial lifecycle patterns were undetected for P. hallaisanensis in the study area. Accordingly, a flowering event of P. hallaisanensis occurred biennially at several microhabitats (other than MH-1 and MH-2), when they contained only first-year seedlings or second-year adults.

The total number of P. hallaisanensis individuals in the study area was 170, 187, and 276 in 2021, 2022, and 2023, respectively ( Figure 2B ). First-year seedlings and second-year adults accounted for 44.4%–87.3% and 12.7%–55.6% of the P. hallaisanensis population, respectively. EYF consistently showed a larger number of P. hallaisanensis individuals than OYF regardless of age throughout the monitoring period ( Figure 2B ).

Sequencing of P. hallaisanensis provided 100.9 Gb of Illumina short reads, which were used to assemble the reference draft genome, involving 2.6 million contigs with a total length of 1.23 Gb and an N50 length of 0.54 Mb, for further SNP identification ( Supplementary Table S1 ). Subsequently, a total of 123.7 Tb of Illumina short reads was obtained from 10 OYF and 10 EYF samples (20 genotypes), from which 3,716 filtered SNPs were identified to analyze the genetic structure across flowering groups and microhabitats. Genetic diversity indices for 3,716 SNPs are described in Supplementary Table S2.

PERMANOVA on the 3,716 SNPs demonstrated that the flowering group (p < 0.005), microhabitat (p < 0.001), and their interaction (p < 0.01) had a significant effect on the multivariate centroid of 3,716 SNPs from P. hallaisanensis ( Figure 3A ). These three factors explained 44.8% of the multivariate variances in SNPs, while the remaining 55.2% of the variance were attributed to within-group variabilities. On the other hand, PERMDISP indicated that only microhabitat (p < 0.05) had a significant effect on the multivariate dispersion of 3,716 SNPs from P. hallaisanensis ( Figure 3A ).

NMDS ordination showed a similar pattern with PERMANOVA and PERMDISP (stress value: 0.12, Figure 3B ). Additional GLM on NMDS axes found that axes 1 and 2 represented the variabilities resulting from microhabitat (p < 0.001) and flowering group (p < 0.005), respectively ( Figure 3B ). Since the multivariate centroids of OYF and EYF were closer in MH-1 than in MH-2 ( Figure 3B ), the multivariate dispersion conversely became larger in MH-2 regardless of the similar genetic diversity indices among the four combinations of OYF and EYF in MH-1 and MH-2 ( Supplementary Table S2 ).

The silhouette width approach suggested that the optimum K number was four, followed by three and two ( Figure 4A ). Neighbor-joining phylogenetic tree and k-means clustering revealed that the 20 genotypes of P. hallaisanensis were divided by microhabitat (MH-1 and MH-2) at K = 2. At K = 3, the genotypes in MH-2 were subdivided by flowering group (OYF and EYF) ( Figure 4B ). Only one of the genotypes of OYF in MH-2 was separated as an additional cluster at K = 4, while the subdivisions between OYF and EYF in MH-1 remained relatively unclear ( Figure 4B ). The STRUCTURE analysis exhibited a similar pattern with the neighbor-joining phylogenetic tree, including the differentiation between MH-1 and MH-2 at K = 2 ( Figure 4C ), and a clearer subdivision by the flowering group in MH-2 than in MH-1 at K = 3 and 4 ( Figure 4C ).