Decoding the genome of Brainea insignis reveals insights into fern evolution and conservation

Domain: biology
Published: 2025-12-30
Editorial outcome: major_revision

Abstract
Ferns are an ancient lineage of vascular plants, yet limited genomic resources constrain both evolutionary and conservation inference. Here, we generate a chromosome-level genome assembly for the endangered cycad fern Brainea insignis (8.62 Gb), the sole species in its genus within eupolypods II, and integrate comparative and population genomics to resolve its evolutionary history and vulnerability. The genome retains the ancient whole-genome duplication shared by leptosporangiate ferns; however, its exceptional size is driven primarily by recent repeat accumulation and further shaped by lineage-specific evolutionary signatures linked to functional specialization. Resequencing across the range identifies three geographically and environmentally structured lineages shaped by Quaternary refugia, limited postglacial expansion and localized admixture. Recently reduced populations show pronounced genomic erosion, including inbreeding and elevated genetic load, due to insufficient time for purging. We detect climate-associated local adaptation and project substantial future genetic offsets, with southwestern Indochina populations at highest risk. Our results expand fern genomics and support spatially tailored conservation strategies that maintain habitat connectivity and promotes adaptive gene flow.

Introduction
Ferns are among the earliest lineages of vascular plants1 and encompass over 13,000 species distributed across diverse ecological niches worldwide2. These species exhibit striking morphological and ecological variation—from small, ground-hugging plants to towering tree ferns—and thrive in habitats ranging from tropical rainforests to arid or temperate zones3,4. As the sister lineage to seed plants, ferns provide critical insights into macroevolution transitions in plant anatomy, life history, and reproduction1,5. Moreover, their high spore production and effective wind dispersal6 make many fern species ideal models for investigating gene flow, dispersal patterns, and local adaptation across broad spatial and ecological gradients.

Despite significant progress in seed plant genomics, which has illuminated the complex genomes of many economically and ecologically important species7–9, fern genomics remains comparatively understudied. This disparity is largely attributable to the large genome sizes and complex ploidy levels characteristic of many ferns10,11. For instance, ferns commonly possess large genomes and intricate chromosome structures, as exemplified by Tmesipteris oblanceolata, which has a record-breaking 160-Gb genome12. Additionally, most population genetic theories and analytical tools are tailored to diploid species13,14, limiting their applicability to ferns, which often deviate from these genetic and genomic norms10,15. As a result, fewer than ten fern species currently have fully assembled genomes, leaving large gaps in our understanding of fern diversification, adaptation, and evolutionary trajectories.

Within ferns, eupolypods—subdivided into eupolypods I (Polypodiineae) and eupolypods II (Aspleniineae)—represent roughly two-thirds of extant fern diversity16. Despite their ecological and evolutionary significance, no high-quality reference genome has been established for this clade. Brainea insignis, commonly known as the cycad fern, belongs to eupolypods II and is the sole species in its genus. This monotypic genus is found in tropical Asia and bears considerable ornamental and medicinal importance17–19. It has been listed as a nationally protected species in China since 199920 and identified as a priority protected species in India21. In China, B. insignis is highly sensitive to environmental disruption and human activities, which have led to precipitous population declines17. Given the species’ dual importance to conservation and evolutionary studies, understanding its genomic makeup and evolutionary resilience is critical for guiding the conservation of both the species itself and the biodiverse habitats it sustains22.

In this study, we present chromosome-level genome assembly for B. insignis, providing an essential genomic resource for the eupolypods II clade. By resequencing 94 individuals from multiple populations, we investigate patterns of genetic diversity, population structure, local adaptation, and the demographic history of this endangered species. Building upon these genomic data, we further estimate genetic offset under future climate scenarios to assess the risks of climate-induced maladaptation. Our results bridge a critical gap in fern genomics, offering fresh insights into the mechanisms underpinning fern genome evolution and informing evidence-based conservation strategies. By addressing key questions regarding the evolutionary potential of B. insignis, our work not only contributes to the broader understanding of fern diversification but also highlights the importance of genomics-based strategies for preserving biodiversity in a rapidly changing world.

Results
Chromosome-scale genome assembly and annotation
To guide our sequencing strategy, we first determined the genome size and chromosome complement of B. insignis. Four independent flow cytometry analyses estimated the genome size at 8.71 Gb (Supplementary Fig. 1), while cytological examination revealed a somatic chromosome number of 2n = 68 (Supplementary Fig. 2). A subsequent genome survey suggested a genome size of approximately 8.40 Gb, with low heterozygosity (0.28%), a high repeat content (89.23%), and diploidy, as indicated by k-mer frequency analyses (Supplementary Fig. 3). Building on these findings, we generated a deep-sequencing dataset to account for the large genome size and the abundance of repeats. We obtained a total of 359.99 Gb (41.76× coverage) of PacBio HiFi reads and 1404.44 Gb (162.92× coverage) of Hi-C reads (Supplementary Table 1 and 2). We then assembled an 8.62 Gb genome, with a contig N50 of 4.36 Mb and a scaffold N50 of 265.61 Mb (Table 1). We anchored 99.81% of the contig length to 34 pseudochromosomes (Fig. 1, Table 1, and Supplementary Fig. 2), corresponding to the 34 chromosomes of the B. insignis haploid set (n = 34). This represents one of the largest haploid genomes with a chromosome-level assembly reported for a non-seed plant, surpassing many existing seed plant assemblies in size.

We assessed the quality of the assembled genome using multiple approaches. First, 99.66% of Illumina short reads (excluding supplemental and secondary reads) mapped to the assembly. Second, BUSCO analysis (viridiplantae_odb12 dataset; updated July 1, 2025) indicated that 97.4% of the 822 conserved genes were completely recovered (Supplementary Table 3). Third, the assembly achieved a LTR Assembly Index (LAI) score of 10.34, surpassing the reference-level threshold of 10. Finally, Clipping Information for Revealing Assembly Quality (CRAQ) analysis showed strong structural accuracy, with a regional assembly quality (R-AQI) of 96.17 and a structural assembly quality (S-AQI) of 98.73. These results collectively underscore the high contiguity, consistency, and completeness of our B. insignis assembly.

Repetitive elements, particularly transposable elements (TEs), account for a large fraction of fern genomes23,24, and B. insignis is no exception. Overall, 7.06 Gb (81.77%) of the genome comprises repeat sequences (Supplementary Table 4), aligning with the genome survey estimate of 89.23% repetition (Supplementary Fig. 3). Long terminal repeat retrotransposons (LTR-RTs) dominate this repeat landscape, accounting for 47.99% of the genome, followed by DNA transposons (DNATs) at 24.49%. The remainder includes non-LTR retrotransposons (e.g., LINEs and SINEs) and a small fraction (2.8%) of unclassified repeats (Supplementary Table 4).

We predicted a total of 43,573 protein-coding genes within the assembled genome (Supplementary Fig. 4a; Supplementary Table 5), with 89.78% of these genes having functional annotations in six major databases (Supplementary Fig. 4b and Supplementary Table 6). The structural characteristics of these genes are broadly consistent with those reported for other fern species (Supplementary Fig. 5): on average, each gene is 15,096.24 bp in length, contains 3.84 exons with a mean exon length of 282.76 bp, and has introns spanning 4,939.39 bp (Supplementary Table 5). Additional BUSCO analysis showed that 692 (84.2% of 822) viridiplantae_odb12 genes were present as complete genes in the annotation (Supplementary Table 3). Collectively, these metrics illustrate the robust quality of our structural and functional gene predictions for B. insignis.

Genomic structural features and comparative genomics
We first examined the structural attributes of the B. insignis genome by profiling gene density, GC content, and the distribution of transposable elements (TEs), particularly long terminal repeat retrotransposons (LTR-RTs) from the Gypsy and Copia families (Fig. 1). Both genes and repeats showed a relatively uniform distribution across the genome, in contrast to the more localized patterns typically observed in seed plants (e.g., Glycine max and Miscanthus floridulus), where gene density generally increases near chromosome termini25,26. Similar homogeneity in genomic architecture has been reported for other homosporous fern genomes, including Ceratopteris richardii27, Alsophila spinulosa28, and Adiantum capillus-veneris29, as well as in select lycophytes (Isoetes taiwanensis30, Huperzia asiatica31, Diphasiastrum complanatum31). Although ferns and seed plants share a common ancestor, these findings suggest that lycophytes and ferns may exhibit more similar genomic structural patterns than seed plants. Accordingly, the differences in gene density and repeat coverage between seed-free and seed-bearing vascular plants may be more nuanced than previously assumed24,29.

Whole-genome duplication (WGD) is a pivotal force shaping genome architecture and driving evolutionary innovation15. Synonymous substitution rates (Ks) distributions are constructed by calculating the Ks between pairs of homologous genes. To explore the role of WGD in B. insignis, we employed a rate-corrected Ks distribution32, which allowed the rescaled orthologous divergence times to be comparable with the paranome Ks distribution of the focal species. A prominent Ks peak at 1.77 suggests an ancient WGD event predating the divergence of Ceratopteris richardii, Adiantum capillus-veneris, Alsophila spinulosa, and Marsilea vestita (Fig. 2a). Further analysis of paralogous Ks distributions and collinearity corroborated this ancient WGD33, pinpointing a consistent peak (~1.7) and revealing no subsequent, lineage-specific WGDs (Supplementary Figs. 6 and 7). These findings align with genome analyses of A. capillus-veneris, which also exhibit only this ancient WGD shared among core leptosporangiate ferns29.

Interestingly, the Ks peak value for B. insignis (~1.7) is lower than that of A. capillus-veneris (~2.15), implying a comparatively slower synonymous substitution rate in B. insignis. To investigate this further, we quantified synonymous (dS) and nonsynonymous (dN) substitution rates in 16 fern species. B. insignis and its close relatives (Woodwardia prolifera) showed significantly lower dS and dN values than many other ferns (Supplementary Fig. 8). Indeed, relative rate tests indicate that B. insignis exhibits a significantly slower evolutionary rate compared to most other core leptosporangiate ferns (Supplementary Table 7), even slower than certain tree fern species (e.g., Alsophila spinulosa). Whole-genome collinearity with the tree fern A. spinulosa revealed extensive one-to-one syntenic blocks in B. insignis (Supplementary Fig. 9). Despite over 200 million years of divergence (TimeTree 5), the breadth of this conserved collinearity indicates that the B. insignis genome has evolved comparatively slowly, at least in terms of genomic architecture. Although B. insignis exhibits a higher ω (dN/dS) ratio than most core leptosporangiate ferns (Supplementary Fig. 10), this may reflect relaxed selective pressures—possibly owing to a stable habitat or small effective population size—rather than accelerated protein evolution. In turn, this slow evolutionary rate may have bolstered genomic stability, enabling B. insignis to persist as the sole species in its genus, although it could also limit its adaptive potential in the face of rapid environmental change.

In addition to WGD, variations in LTR-RT content strongly influence genome size evolution in many plant lineages24,34,35. However, their specific impact on fern genome expansion remains poorly understood, largely due to limited genomic data. To address this gap, we compared LTR-RT composition and insertion times in eight fern genomes. Our findings indicate that larger fern genomes tend to harbor a higher proportion of LTR-RTs with earlier insertion time (Fig. 2b; Supplementary Fig. 11). While these observations suggest a positive correlation between LTR-RT activity and genome size expansion in ferns, the limited number of available fern genomes hampers broad statistical inferences. Consequently, more extensive sampling and high-quality assemblies are needed to better understand how LTR-RT dynamics have shaped the remarkably diverse genome architectures across ferns.

To further investigate gene family evolution in the B. insignis genome, we constructed a high-confidence phylogeny for 16 species using 103 single-copy gene families (Fig. 2c). We identified 353 significantly expanded and 200 significantly contracted gene families on the branch leading to B. insignis, with the expanded families enriched in biological processes related to cell wall organization and lignin metabolism (“plant-type cell wall organization”, “lignin metabolic process” and “lignin catabolic process”; Fig. 2d). Gene families involved in monolignol biosynthetic pathway exhibit Ka/Ks < 1 in B. insignis and two other tree ferns (Supplementary Fig. 12), and microsynteny analyses further reveal that a subset of these lignin-related genes is highly conserved across species (Supplementary Fig. 13). These findings suggest that the ancestral lineage of B. insignis evolved specialized lignin-related traits—such as robust, lignified structures—that have been retained in the modern species, while key metabolic pathways remained long-term functional stability. Interestingly, a recent study showed that Stenochlaena palustris, a member of the same family, possesses remarkable lignin architectures36, indicating that S-lignin production evolved independently in ferns. Although lignin composition in B. insignis remains to be characterized, the observed expansion and synteny suggest lineage-specific elaboration of lignin-related capacities. Together, these observations are consistent with lineage-specific diversification of lignin pathways in ferns.

Overall, our findings demonstrated that an ancient WGD, repetitive-element dynamics, and gene family expansions have collectively sculpted the genome of B. insignis. Although B. insignis shares a deep polyploidization event with other core leptosporangiate ferns, it exhibits a slower evolutionary rate and a larger genome, underscoring the complexity of diploidization processes and the importance of repetitive elements in shaping fern genome diversity.

Population structure and demographic history
We sequenced 94 B. insignis individuals from 29 geographic locations to an averagedepth of 21.81×, generating 17.79 Tb of raw data (Fig. 3a; and Supplementary Fig. 14 and Supplementary Data 1). After mapping paired-end reads to the B. insignis reference genome and applying stringent filtering criteria, we identified 75,060,153 high-quality SNPs and 9,414,109 core variants defined as a high-confidence SNPs set retained after stringent quality control (Supplementary Fig. 15). Dataset 1-3 were derived from these variants and showed no evidence of chromosomal bias (Supplementary Fig. 16).

To elucidate the genetic structure of B. insignis, we applied several methods, including ADMIXTURE, principal components analysis (PCA), neighbor-joining (NJ) trees, and chloroplast haplotype networks. ADMIXTURE identified two primary clusters (K=2), differentiating Yunnan populations (YN lineage) from those in southern China (SC lineage). The optimal number of clusters was determined to be three (K=3), revealing a distinct lineage from Vietnam (VN lineage) and two apparent admixture zones (Admixture1 and Admixture2) (Fig. 3b; and Supplementary Fig. 17). Additional analyses corroborated this tri-lineage pattern: (1) the chloroplast haplotype network revealed clear distinctions among the YN, VN, and SC lineages, with apparent admixture haplotypes (Supplementary Fig. 18); and (2) PCA and NJ trees confirmed genetic separations corresponding to these three main lineages while highlighting intermediate signatures in the admixture zones (Supplementary Fig. 19). Together, these results show strong isolation by distance (IBD) and isolation by environment (IBE), suggesting that both geographic distance and environmental heterogeneity jointly influence the spatial distribution of genetic variation. This pattern is consistent with multiple glacial refugia during the Quaternary, followed by post-glacial expansions and secondary contact among lineages.

We reconstructed the demographic trajectories of each lineage using the pairwise sequentially Markovian coalescent (PSMC) method (Fig. 3c). All lineages exhibited similar initial trends: an increase in effective population size (Ne) peaking around 2 million years ago (Mya), followed by a protracted decline until approximately 300 thousand years ago (Kya). A subsequent expansion culminated in a sharp bottleneck around 200 Kya, coinciding with the onset of Marine Isotope Stage 6 (MIS6) glaciation37,38. These observations mirror the glacial-contraction and interglacial-expansion cycles reported for other East Asian relict species, including Ginkgo biloba39, Cercidiphyllum japonicum40, and Davidia involucrata41. After the MIS6 bottleneck, the YN lineage diverged from the common ancestor around 100 Kya, whereas the VN and SC lineages separated around 40 Kya. Each lineage subsequently maintained relatively small and stable Ne values up to the present (Fig. 3c; and Supplementary Fig. 20).

To complement the PSMC analyses, we performed coalescent-based simulations for each lineage using fastsimcoal2 and identified the most likely speciation model among five initial scenarios without introgression (Supplementary Fig. 21 and 22; Supplementary Table 8 and 9). After determining the best-fitting model, we tested nine additional models incorporating gene flow (Supplementary Fig. 23; and Supplementary Table 10). The optimal model corroborated the PSMC-derived divergence times (Fig. 3d) and provided a robust framework for understanding the evolutionary history of B. insignis, further supporting its status as a relict fern species17,20. Our simulations also revealed variations in gene flow across different divergence phases. Substantial sharing of identical-by-descent haplotypes between the YN/SC lineages and the admixture zones (Admixture1 and Admixture2) implies relatively recent inter-lineage gene flow (Fig. 3e). Furthermore, TreeMix analysis, with an optimal migration edge of 1, indicates recent gene flow from YN to VN (Supplementary Fig. 24), underscoring ongoing genetic exchange and secondary contact among geographically proximate populations.

Genomic effects of recent population decline on inbreeding and genetic load
Although effects of population declines are well-documented in many seed plants and animal lineages, the genomic consequences of population decline in ferns remain largely unexplored42. To address this gap, we analyzed genetic diversity, inbreeding, and genetic load in B. insignis. We hypothesized that the consistently low effective population size (Ne) since the last bottleneck, compounded by severe declines in recent decades, has led to further losses of genetic diversity. Indeed, genome-wide heterozygosity in the three lineages (YN, VN, and SC) remains low, averaging 0.08, 0.12, and 0.10, respectively (Fig. 4a). Accordingly, the YN lineage exhibits the lowest nucleotide diversity (mean π = 1.043×10−3), followed by VN (1.126×10−3) and SC (1.379×10−3) (Supplementary Fig. 25). YN also exhibits slower linkage disequilibrium (LD) decay, which is consistent with its lowest level of genetic diversity. Genetic divergence (FST) values (0.116–0.320) correlate with the geographic distances between the lineages, while positive Tajima’s D values corroborate the demographic bottlenecks inferred from our historical reconstructions.

To quantify inbreeding, we identified runs of homozygosity (ROH) in each lineage. All three lineages show numerous long ROHs (>100 kb), with the fraction of the genome in ROH (FROH) averaging 0.51 (YN), 0.50 (VN), and 0.43 (SC) (Fig. 4b). Furthermore, a substantial proportion of ROHs exceed 1 Mb (Fig. 4c), indicating pronounced inbreeding. We observed a strong negative correlation between individual heterozygosity and FROH (R2=0.77, P<0.001), confirming that high levels of inbreeding are reducing genome-wide diversity (Fig. 4d).

Next, we assessed genetic load by examining allelic states at loss-of-function (LOF) and missense (deleterious) variants. The ratio of nonsynonymous to synonymous (π0/π4) diversity was relatively high, ranging from 0.431 to 0.503 (Supplementary Table 11), suggesting that purifying selection on 0-fold sites has been weak or relaxed in B. insignis. To put our estimates in a broader context, we conducted a comparative analysis of π₀/π₄ ratios across more than 40 plant species, spanning ferns and seed plants (Supplementary Fig. 26). The endangered ferns (B. insignis and Alsophila spp.) exhibit markedly elevated π₀/π₄ ratios relative to most seed plants, with the notable exception of Acer yangbiense, a species characterized by an extremely small population size. These results support the expectation that small population sizes weaken purifying selection and inflate genetic load across disparate plant lineages. Although absolute values may vary with sampling and analytical pipelines, the qualitative pattern is robust.

Deleterious (DEL) and LOF mutations in the homozygous state—an indicator of genetic load43,44—are most frequent in the YN lineage (Fig. 4e, f). Moreover, applying a Grantham Score threshold of 150 further underscores the elevated fraction of deleterious missense variants in YN lineage (Supplementary Fig. 27). The frequency of homozygous DEL and LOF alleles correlates with FROH (Fig. 4g, h). We also observed a higher density of LOF variants within ROH segments than outside them (Fig. 4i), highlighting the role of inbreeding in exacerbating genetic load. Functional annotations link LOF-afflicted genes to stress resistance and DNA repair (Supplementary Data 2)—processes vital for organismal survival—implying that accumulated deleterious mutations may diminish population fitness.

Collectively, our findings demonstrate that B. insignis harbors extremely low genetic diversity in all three lineages and exhibits pronounced inbreeding and genetic load, particularly in the YN lineage. This pattern contrasts with long-term small populations (e.g., Alsophila spinulosa)42, where signals of purging have been reported43,45–48. In B. insignis, however, a recent and severe population contraction likely limited purging opportunities, thereby facilitating the accumulation of deleterious alleles. Consequently, these results underscore that fern lineages shaped by different demographic histories can display sharply contrasting patterns of genetic load purging. In the case of B. insignis, the elevated genetic load presents serious challenges to its long-term viability, emphasizing the need for conservation measures that safeguard both habitat protection and the genetic health of the remaining populations.

Adaptive differentiation and genetic vulnerability under future climate change
To identify genes potentially involved in lineage-specific adaptations, we conducted selective sweep analyses on the southern China (SC) and Yunnan (YN) lineages, which had large sample sizes and exhibited the highest genetic differentiation. We detected 3,846 positively selected sites (involving 116 genes) in the SC lineage and 2,299 such sites (encompassing 64 genes) in the YN lineage (Fig. 5a; and Supplementary Fig. 28). Only three of these genes were shared between the two lineages. Gene Ontology (GO) analysis revealed that most positively selected genes were linked to core cellular and metabolic processes, including binding and catalytic functions (Fig. 5b). Although the SC and YN lineages shared relatively few positively selected genes, they exhibited considerable overlap in functional categories, suggesting functional convergence on likely adaptive traits. Nonetheless, differences in gene counts per category and lineage-specific annotations point to diverging selection pressures that have driven adaptive differentiation.

Next, we examined the relationship between allele frequency and environmental variables in the context of local adaptation. We identified 4,974 core adaptive variants (Fig. 5c) that displayed significant isolation-by-distance (Mantel’s r=0.4301, P<0.05) and isolation-by-environment (Mantel’s r=0.4058, P<0.05) patterns (Supplementary Fig. 19). However, after controlling for geographic and environmental effects separately, neither relationship remained significant (partial Mantel tests: P=0.147 for IBE, P=0.179 for IBD), suggesting that both factors jointly shape adaptive genetic variation. Redundancy analysis (RDA) identified three distinct genetic clusters corresponding to the YN, SC, and VN lineages, with admixture populations occupying intermediate positions (Supplementary Fig. 29). The YN lineage correlated with a high diurnal temperature range (bio2), while the SC lineage was linked to low bio2, high annual mean temperature (bio1), and high precipitation in the driest month (bio14). In contrast, the VN lineage formed its own distinct environmental cluster, and the admixture populations displayed heterogeneous environmental signatures. Collectively, these findings support local adaptation in B. insignis, driven by the interplay of environmental variation and geographic distance.

To evaluate the risk of genetic maladaptation under global warming, we employed a genetic offset framework49,50. Under various emission scenarios, higher greenhouse gas concentrations consistently produced larger offsets, signaling a greater risk of maladaptation (Supplementary Fig. 30). Populations in the southwestern Indochinese Peninsula exhibited relatively high local offsets, whereas populations outside this region showed lower offsets (Supplementary Fig. 31a, d). Forward and reverse offsets were generally very low across the species’ range, suggesting that future habitats may remain suitable (Supplementary Fig. 31b, c, e, f). However, these theoretical predictions must be interpreted in light of dispersal limitations, as B. insignis spores typically travel only tens to hundreds of kilometers51. Under high genetic offsets (>0.15), most populations are strongly locally adapted to future climates, with only a few in the Indochinese Peninsula requiring long-distance migration (Supplementary Fig. 32a). Considering moderate offsets (>0.05), the required migration distances increase by several thousand kilometers, with many populations needing to disperse over 5,000 km to avoid maladaptation, far exceeding their natural dispersal capacity (Supplementary Fig. 32b).

Collectively, these results indicate that populations in the southwestern Indochinese Peninsula face elevated genetic vulnerability and high risk under future climate change (Fig. 5d). Polar plotting revealed no clear directional trend for potential migration (Supplementary Fig. 33), reinforcing the likelihood that range fragmentation, rather than large-scale migration, will shape the future distribution of B. insignis. Overall, these findings highlight not only the magnitude of genetic vulnerability but also the limited evolutionary and dispersal capacity of B. insignis under rapid climate change. The pronounced regional differences in offset levels suggest that conservation strategies may need to be spatially tailored, prioritizing populations in the southwestern Indochinese Peninsula. Moreover, the extreme migration distances implied by moderate offsets emphasize that natural dispersal alone is unlikely to maintain adaptive potential, raising the possibility that assisted gene flow or habitat connectivity restoration may be necessary to safeguard long-term persistence.

Methods
Plant materials and genome sequencing
Fresh fronds from a single B. insignis individual were collected from the South China National Botanical Garden in Guangzhou, China. Total genomic DNA was extracted using the CTAB method52. A preliminary genome survey was performed on the BGI DNBSEQTM platform, generating 512.37 Gb of short-read data to estimate genome size, heterozygosity, and repeat content. High-molecular-weight DNA libraries (15–18 kb) were constructed and sequenced on the PacBio Sequel II/IIe platform in CCS (HiFi) mode, generating 359.99 Gb of HiFi reads with an average read length of 16.54 kb. Hi-C libraries were prepared according to established protocols53,54. Briefly, fresh fronds were fixed with 4% formaldehyde, digested with DpnII, and biotin-labeled DNA ligation products were sheared and used to build Illumina PE150 libraries. Five Hi-C libraries were sequenced; each amplified for 12–14 PCR cycles.

Genome size estimation and chromosome counting
A 17-mer frequency distribution was generated using KmerFreq (v4.0 in GCE v.1.0.2), and the genome size was calculated in GCE (v.1.0.2)55. Flow cytometry was performed for confirmation: leaf nuclei from B. insignis and the internal control Camellia sinensis var. assamica56 were co-stained with propidium iodide and analyzed using a BD FACScalibur flow cytometer at 488 nm excitation. The C-values were calculated based on the PI fluorescence peaks of the sample and control. For chromosome counting, root tips (~1 cm) from the sequenced individual were pretreated with p-dichlorobenzene for 5 h at room temperature, fixed in 3:1 (v/v) ethanol-to-acetic-acid solution at 4 °C, then macerated in 1 M HCl at 60 °C for 8 min. Chromosomes were stained with carbolfuchsin and observed by optical microscopy57.

Genome assembly and quality assessment
Genome assembly was performed with Hifiasm (v.0.19.5-r587)58 in Hi-C mode using default parameters, incorporating both HiFi and Hi-C reads. The primary assembly was 8.78 Gb in length, consisting of 5,791 contigs (contig N50 of 4.36 Mb; Table 1). Hi-C reads were mapped to the primary assembly with BWA (v.0.7.17-r1188)59, scaffolding was performed using Juicer (v.1.6)60 and 3D-DNA pipeline61, followed by manual curation using Juicerbox module.

Assembly quality was evaluated by mapping short reads with BWA-MEM, completeness was assessed with BUSCO (viridiplantae_odb12 dataset)62. Additionally, N50 and other quality metrics were obtained using QUAST (v.5.2.0)63. Finally, CRAQ (v.1.0.9)64 was used to identify regional (R-AQI) and structural (S-AQI) assembly errors.

Genome annotation
We employed an integrated approach of homology alignment and de novo search to annotate repetitive elements. RepeatMasker (http://www.repeatmasker.org) with Dfam65 and Repbase66 databases was employed for homology-based detection. For ab initio prediction, a de novo transposable element (TE) library was generated using MITE-Hunter67 and RepeatModeler (v.2.0.3)68. The unknown elements were classified by DeepTE69, and a non-redundant TE library was produced via Uclust70. Final repeats were masked with RepeatMasker.

Protein-coding genes were predicted using ab initio, homology-based, and RNA-seq–assisted approaches (details in Supplementary Note 1). Functional annotations were assigned via BLASTp (E-value ≤ 1e−5) against Swiss-Prot71, with motifs and domains identified by InterProScan (v.5.39)72. Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway annotations were assigned accordingly. Non-coding RNAs were characterized as described in the Supplementary Note 2.

Comparative genomics analysis
Synteny and whole-genome duplication events were evaluated using WGDI (v0.6.1)33 for collinearity and Ks distributions. Substitution rate variation was accounted for by employing four species trios: (1) B. insignis, Adiantum capillus-veneris, Lygodium japonicum; (2) B. insignis, Ceratopteris richardii, L. japonicum; (3) B. insignis, Alsophila spinulosa, L. japonicum; and (4) B. insignis, Marsilea vestita, L. japonicum. We standardized the synonymous substitution rates for each trio from divergence events between species to the focal species (B. insignis)73.

For analyses of LTR retrotransposons (LTR-RTs), we obtained genomic data from seven available fern species—Salvinia cucullata74, Azolla filiculoides74, Marsilea vestita75, Adiantum capillus-veneris29, Alsophila spinulosa28, Adiantum nelumbodies76, and Ceratopteris richardii27—from publicly available datasets. LTR-RTs were identified uniformly in eight fern species using LTR_Finder77, LTRharvest78, and LTR_retriever79. We applied a phylogenetically informed calibration to account for lineage-specific rate heterogeneity80. For LTR retrotransposons, we estimated the corrected nucleotide substitutions (K) between paired 5’ and 3’ LTRs using LTR_retriever, which employs the Jukes-Cantor (JC69) model. The resulting K values were converted into LTR insertion times using lineage-specific substitution rates calibrated on a phylogenetic framework, enabling cross-species comparisons. The results were visualized in a density plot using ggplot2 (v.3.5.1)81.

To investigate evolutionary rates, we performed phylogenomic analyses on 16 representative fern species using 8,720 low-copy and 31 single-copy orthologs identified by Orthofinder (v.2.5.5)82. Multiple sequence alignments were performed in MAFFT (v.7.520)83, trimmed by trimAL (v.1.4.rev15)84. A maximum likelihood (ML) phylogeny was then inferred with IQ-TREE (v.2.2.6)85, employing 1,000 bootstrap replicates. Substitution rates (dS, dN) were estimated with PAML (v.4.8)86, and a likelihood ratio test compared free-ratio and one-ratio models. The significance of differences in synonymous (dS) and nonsynonymous (dN) substitution rates among the predefined groups (single-copy and low-copy) was evaluated using analysis of molecular variance (AMOVA). Finally, relative rate tests (RRT) were performed in MEGA87 on a concatenated protein alignment of 31 single-copy genes.

Gene family identification and evolution
For gene family evolution, a phylogenetic tree was constructed across 16 taxa, including a bryophyte (Physcomitrium patens), a lycophyte (Selaginella moellendorffii), four seed plants, and ten ferns (including B. insignis). We rooted the resulting tree with P. patens, converted it into an ultrametric chronogram using r8s88, and calibrated it against TimeTree (v.5.0)89. CAFE (v.4.2.1)90 was used to identify significantly expanded or contracted gene families (P<0.05). We characterized 11 gene families associated with the monolignol biosynthesis based on homology to Alsophila spinulosa genes. One-to-one orthologs in B. insignis, A. spinulosa, and Sphaeropteris lepifera were identified by reciprocal BLAST, and Ka/Ks ratios were calculated in KaKs_Calculator (v.3.0)91 to assess selective pressures on lignin-related genes.

Resequencing, SNP calling, and filtering
Fresh leaves from 94 B. insignis individuals representing 29 geographic locations (Supplementary Data 1) were used for DNA extraction using the CTAB method52. Paired-end libraries were sequenced on the DNBSEQ-T7 platform, and raw reads were quality-filtered with fastp (v.0.22.0)92. Clean reads were aligned to the B. insignis reference genome using BWA-MEM59. The resulting SAM files were converted to BAM format and sorted with SAMtools93. PCR duplicates were removed via Picard (https://broadinstitute.github.io/picard/). Variants were called following the GATK (v.4.5.0.0) pipeline94. Low-confidence variants were removed using the following strict hard-filtering thresholds: QD<2.0 || MQ<40.0 || FS>60.0 || SOR>3.0 || MQRankSum<-12.5 || ReadPosRankSum<-8.0. Subsequently, a multi-step soft-filtering process was applied to generate multiple high-quality SNP datasets for downstream analyses (Supplementary Fig. 15). We define the high-confidence variants (core variants) as biallelic SNPs retained after joint genotyping and stringent hard/soft filtering. All downstream variant panels (Dataset 1-3) were derived from this core variant set (Supplementary Fig. 15).

Population structure and genetic diversity analyses
To assess population structure, we generated 342,582 unlinked SNPs (minor allele frequency > 5%), which were analyzed with ADMIXTURE (v.1.3.0)95 for K=1 to K=8 for cross-validation. Principal components analysis (PCA) was performed in PLINK (v.1.9)96, and a neighbor-joining (NJ) phylogenetic tree was constructed using PHYLIP. A haplotype network was inferred with PopART97 using the 95% statistical parsimony (TCS) network method (Supplementary Note 3). Nucleotide diversity (π), genetic differentiation coefficient (FST) and Tajima’s D were estimated in 100-kb windows using pixy (v.1.2.10)98 and VCFTools (v.0.1.16)99. Linkage disequilibrium (LD) decay was calculated via PopLDdecay (v.3.42)100.

Demographic history inference
For demographic history, we selected six samples (≥20× coverage) per major lineage to build diploid consensus sequences using BCFtools (v.1.14)101. PSMC102 was run with default parameters (25 iterations, -N25, -r5, -p 4+25*2+4+6) under a mutation rate of 1.25×10−8 per site per generation and a 5-year generation time. In addition, SMC++ (v.1.15.2)103 was performed to infer a finer-scale reconstruction of recent demographic histories within each lineage. To model complex demographic scenarios, we employed a stepwise, hierarchical approach using fastsimcoal2104. We initially modeled simpler histories with fewer parameters (e.g., divergence without migration), and progressively incorporated additional features such as continuous migration and pulse admixture. A total of 10 demographic and 5 differentiation models (Supplementary Fig. 21 and 22) were evaluated, each with 100 independent runs. The best-fitting scenario was selected based on the lowest AIC and ΔLhood (Supplementary Table 8–10). This stepwise strategy balances model realism with statistical robustness and mitigates potential identifiability issues arising from simultaneous estimation of many parameters.

To investigate gene flow, we generated a maximum likelihood drift tree using TreeMix (v.1.13)105 with migration edges from 1 to 5, selecting the optimal number with OptM106. Finally, identity-by-descent (IBD) segments were identified using the BEAGLE algorithm107.

Estimation of inbreeding and genetic load
Individual heterozygosity was estimated with ANGSD (v. 0.937)108 under a folded site frequency spectrum (fSFS) model. Runs of homozygosity (ROHs) were identified with PLINK (v.1.9)96 using a minimum of 50 consecutive SNPs per 100-kb window. The inbreeding coefficient (FROH) was calculated as the proportion of the genome covered by ROHs109,110.

Derived alleles were defined as those where more than 50% of individuals carried the same homozygous genotype111. SnpEff (v.4.3t)112 was used to annotate missense and loss-of-function (LOF) variants. Missense mutations were further assessed using SIFT-4G113 to classify potentially deleterious (DEL) substitutions, and the Grantham Score (GS≥150)114 was applied to confirm deleterious missense variants. We calculated the proportion of homozygous DEL and LOF sites relative to the total number of homozygous plus heterozygous sites of each individual. Correlations among heterozygosity, inbreeding, and deleterious load (DEL and LOF) were visualized with SRplot115. We compared homozygous LOF variant frequencies inside and outside ROHs to evaluate the impact of inbreeding on genetic load. To quantify nonsynonymous to synonymous diversity, we identified 0- and 4-fold degenerate sites from the reference CDS using the degenotate pipline (https://github.com/harvardinformatics/degenotate), extracted the corresponding population polymorphism data from variant call format (VCF) file, and calculated nucleotide diversity (π) at these sites with ANGSD108.

Selective sweep analysis
Selective sweeps were examined for two major lineages (YN and SC). We used a combined approach incorporating high FST signals (top 5% in 500 kb windows, 5 kb step) and within-lineage sweeps identified by RaiSD (v.2.9)116. Candidate selective sweeps were defined where both methods overlapped, with the 99.99% quantile of μ serving as the threshold. Manhattan plots were generated in R using the qqman package117.

Genotype-environment association analyses
We retrieved 19 bioclimatic layers (1970–2000) from WorldClim (https://worldclim.org/) using ArcGIS v.10.8 for spatial extraction and retained seven uncorrelated variables (|r|<0.7) via Pearson correlation: bio1, bio2, bio3, bio7, bio12, bio14, and bio15 (Supplementary Fig. 34). To identify a robust set of climate-adapted loci, we applied Latent Factor Mixed Models (LFMM) and Redundancy Analysis (RDA) to test genotype-environment associations (GEA). LFMM were conducted using three latent factors in the LEA package (v.3.16.0)118 to account for population structure. Variants significantly associated (FDR<0.05) with at least three environmental variables were considered outliers. RDA was performed using vegan (v.2.6-8)119. Significant variants were identified based on extreme loadings (±3.5 SD, equivalent to a two-tailed P value of 0.0005) along one or more RDA axes.

To disentangle the effects of isolation by distance (IBD) from isolation by environment (IBE), we conducted Mantel and partial Mantel tests separately on adaptive and neutral SNP datasets. Pairwise matrices of FST/(1−FST) were tested for correlation with geographic and environmental distances using the vegan package.

Genetic offset under future climate scenarios
Genetic offset (also termed genomic vulnerability or risk of non-adaptedness, RONA49,50) was estimated specifically using the adaptive SNPs. We applied generalized dissimilarity modeling (GDM)120,121 with future climate projections (2081–2100) from WorldClim 2.1 under two climate models (BCC-CSM2-MR and GISS-E2-1-G) and two Shared Socioeconomic Pathways (SSP245 and SSP585) at 2.5 arc minutes resolution. Because the two climate models showed a strong correlation (Supplementary Fig. 35), BCC-CSM2-MR was used. Forward and reverse genetic offsets120,122 were calculated to account for the potential effects of migration, then normalized across each channel (RGB) for spatial visualization. The local, forward, and reverse offsets were then mapped to the red, green, and blue channels, respectively. For each grid, we calculated genetic offsets and the maximum allowable migration distance and used Matplotlib (https://matplotlib.org/) to create a polar plot representing migration distance and direction.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Xia et al presented the first fern genome from the eupolypod II clade, and characterized the history of whole genome\nduplication, repeat expansion, gene family evolution, and population demography. I found the results on WGD and repeats\nvery interesting, which have broad implications in understanding the gigantic genomes in ferns. On the other hand, I’m\nsomewhat skeptical about authors’ take on the evolution of lignin genes. Without knowing exactly what kinds of unique\nlignin Brainea has, I feel the authors are overinterpreting their CAFÉ results. Furthermore, with their sampling, the\nexpansion/contraction is not specific to Brainea per se, but the >100-million-years-old lineage leading to the species. I’m not\nan expert of population genomic analyses and cannot comment on that part of the manuscript.\nLine-by-line comments below:\nLine 123 and 147. The BUSCO scores are low, but I think it’s because the authors used the embryophyta set, which is\nheavily based on angiosperm genomes. I’d suggest using the viridiplantae set, one that has more representations from non-\nseed plants.\nLine 181-186. I don’t think the data presented support that Brainea has had low substitution rate. From extended fig. 3, I can\nonly see lower rate from Ks on single copy genes, and the rest look all the same. But in the text, it was stated that “B. insignis\nand its close relatives (Woodwardia prolifera and Pronephrium simplex) showed significantly lower dS and dN values than\nother ferns (Extended Data Fig. 3).” I just don’t see the pattern and what stats were done to support this?\nLine 207-223. I think the interpretation here is problematic. The inferred expansions/contractions did not necessarily happen\nto B. insignis specifically, but the lineage leading to it. The closest species on this CAFÉ analysis was >100 million years\napart. Therefore, there is no evidence that the expanded lignin biosynthesis genes “plays a crucial role in the persistence of\nthis monotypic genus”\nLine 330. “To identify genes underlying lineage-specific adaptations” – this to me sounds like you have evidence for local\nadaptations. But I’m not sure if that evidence is presented.\nLine 385-386. I’m not sure I follow how slow substitution reflects “ecological adaptability”.\nFig. 1b. There is a sudden spike in gene density at Chr04? Curious to know what that is about.\nFig. 3. It would be nice to know the distribution range of this species, so that one can assess if the sampling is representative\nand capturing the geography.", "reviewer_2": "To the Authors and Editors,\nXia et al. provide a well-written manuscript of a well-executed study detailing the assembly and annotation of the Brainea\ninsignis genome, as well as a population genomic assessment of the genetic diversity and demographic history of this\nendangered fern species. After several readings of the manuscript and supplemental information, I cannot find any major\nflaws in their study design, methods, or data analyses related to the genome sequencing, assembly, annotation, or\npopulations genomic analyses. The work (with one exception) meets the standards of the field, and I feel that the methods\nare of sufficient detail to be reproducible. I will refrain from assessing their evaluation of genetic offset under future climate\nscenarios–this is too far afield from my expertise. I do, however, take some issue with some of their wording and\ninterpretation of results, and I will detail a few examples at the end of this review.\nThis work has substantial significance for the fields of comparative plant genomics, fern evolution and population genetics,\nand conservation genomics. In my opinion, the greatest impact of this publication is that it provides the first chromosome -\nlevel assembly for the Eupolypod II ferns, representing approximately 20% of extant fern diversity. It fills a major gap in\ngenomic resources for ferns. However, Brainea insignis, a monotypic genus and “living fossil” akin to Gingko biloba and\nCercidiphyllum japonicum, is an odd, atypical choice to represent this group. Perhaps, this choice reflects the authors’ focus\non conservation, rather than comparative genomics for this study?\nNonetheless, this genome assembly is a novel and needed resource for meaningful comparative genomics investigations of\nchromosome number, genome size, genic content, and rates of evolution relative to other ferns and to other major plant\nlineages, especially their sister group the seed plants. The authors briefly touch on this importance throughout the\nmanuscript by confirming an ancestral WGD for the core leptosporangiate ferns, comparing intron sizes/exon numbers\nacross ferns, illustrating the positive relationship between TE content and genome size in ferns, assessing rates of evolution\namong ferns, etc. In my opinion, the appeal and application of this manuscript would be broadened by refocusing and\nexpanding discussion of these topics.\nA second major contribution of this work is a thorough assessment of population genetic structure in a biogeographic\ncontext, demographic history, gene flow, and signatures of local selection. The study of population genomics in ferns,\nparticularly employing whole genome resequencing, is essentially non-existent and has been long hampered by a lack of\ngenomic resources. I can think of only one other study ( Yi et al. 2024. The Plant Journal 120(4):1522-35) that utilized a\ncomparable genomic-scale dataset to estimate signatures of inbreeding and mutational load in a fern. I believe this\nmanuscript will become a crucial resource/protocol for other researchers interested in population-level inquiries in ferns, in\nparticular, for address many long-standing but previously intractable questions related to breeding systems in organisms that\nare spore dispersed and have independent sporophyte/gametophyte life phases.\nA third contribution of this work is that is provides further support for the independent evolution of lignin syringly “S” units,\nmainly reported from angiosperms, in B. insignis. This finding builds on a recent Nature Plants publication (Ali et al. 2024.\nhttps://doi.org/10.1038/s41477-025-684 01978-y), indicating that lignin S-units have evolved repeated across land plants,\nand that ferns have independently evolved mechanisms of structural support and vasculature in terrestrial environments.\nAs mentioned in the first paragraph, there are several instances where I feel the authors inaccurately depicted background\ninformation or have interpreted their results too liberally. Below are three examples, that could be amended:\n1-Line 60 “Due to their phylogenetic position bridging non-vascular bryophytes and seed-bearing vascular plants…” This\nmay seem overly picky, but I think it misrepresents the placement/informativeness of ferns in plant evolution. Ferns do not\n“bridge” or have some type of “intermediate” phylogenetic placement. (As the authors know) ferns are the sister group to\nseed plants, and share a common ancestor with seed plants, not shared with bryophytes. Please clarify for accuracy.\n2-Line 130 “These results collectively underscore the high contiguity, consistency, and completeness of our B. insignis\nassembly.” I am concerned about the authors’ conclusion of completeness given the BUSCO score of 81.5% complete\nrecovery (77.7% based on annotation assessment, Supplementary Table 3). Despite having high LAI and CRAN metrics,\nthis is a very low BUSCO score (The expected standard for the field is usually 95% complete or higher). At minimum, the\nauthors need to address the possible reasons for the low BUSCO recovery and discuss whether it is an artifact of their\nassembly or possibly reflecting a biological reality (unlikely).\n3-Line 187 “relative rate tests suggest that B. insignis has the slowest evolutionary rate among core leptosporangiate ferns\n(Supplementary Table 7).” Given the information in Supplementary Table 7, this is an overexaggeration. It appears that B.\ninsignis does not have a significantly slower rate of evolution (in 31 single copy genes) relative to several other core\nleptosporangiate ferns, including L. chinensis and D. repens. Please amend this statement and temper your conclusion that “\nthis slow evolutionary rate may have bolstered genomic stability, enabling …”.\nThank you for the opportunity to review this manuscript."}, "author_response": "1 Response to Reviewer’s comments\n4 the history of whole genome duplication, repeat expansion, gene family evolution, and\n5 population demography. I found the results on WGD and repeats very interesting, which\n6 have broad implications in understanding the gigantic genomes in ferns. On the other\n7 hand, I’m somewhat skeptical about authors’ take on the evolution of lignin genes.\n9 are overinterpreting their CAFÉ results. Furthermore, with their sampling, the\n10 expansion/contraction is not specific to Brainea per se, but the >100-million-years-old\n11 lineage leading to the species. I’m not an expert in population genomic analyses and\n12 cannot comment on that part of the manuscript.\n14 manuscript. We are glad to hear that you found our results on whole genome duplication\n15 (WGD) and repeat expansion interesting. Regarding the evolution of lignin genes, we\n16 understand your skepticism. We acknowledge that the exact composition of lignin in\n20 lineage leading to B. insignis, rather than changes specific to the species itself.\n22 sampling focused on lineages more closely related to B. insignis (see Response Fig. 1a).\n24 lineage-specific shifts from more ancient trends. Notably, we still detected 353\n25 significantly expanded gene families (p < 0.05) in B. insignis that are enriched in\n26 biological processes related to “cell wall organization” and “lignin metabolic process”\n27 (Response Fig. 1b). In contrast, a parallel analysis in the herbaceous fern Adiantum\n28 capillus-veneris, which lacks lignified woody structures, revealed no enrichment in\n29 lignin-related pathways among its significantly expanded gene families (Response Fig.\n30 1c). These results strengthen our confidence that lignin biosynthesis-related gene\n31 families have indeed undergone expansion in B. insignis, potentially reflecting the\n32 species’ adaptation to a woody, tree-fern growth habit. We have revised the manuscript\n33 to clarify that these expansions are inferred for the lineage leading to B. insignis and to\n34 avoid overstating their importance without direct chemical evidence of lignin\n35 composition in this species (Manuscript Lines 208-211).\n36\n39 the embryophyta set, which is heavily based on angiosperm genomes. I’d suggest using\n40 the viridiplantae set, one that has more representations from non-seed plants.\n42 analysis. As correctly noted, our initial use of the embryophyta_odb10 dataset, which\n43 is biased toward angiosperms and may not be ideal for fern genome. In response, we\n44 have re-evaluated both our genome assembly and annotation using the more inclusive\n45 viridiplantae_odb12 dataset (updated July 1, 2025), which include a broader\n46 representation of non-seed plants. The updated BUSCO results (presented in Response\n48 completeness. For example, the genome assembly’s BUSCO completeness increased\n49 from ~81.5% to 97.4% using the viridiplantae dataset, indicating that the earlier lower\n50 score was indeed an artifact of the reference dataset choice rather than genome quality.\n52 the rigor of our genomic assessment.\n54 viridiplantae_odb12 dataset\nEvaluation Complete Single-copy Duplicated Fragmented Missing Total\nType (C) (S) (D) (F) (M) BUSCOs (n)\nGenome\n97.4% 89.9% 7.5% 2.1% 0.5% 822\nAssembly\nGene\n84.2% 79.1% 5.1% 12.7% 3.2% 822\nAnnotation\n56 substitution rate. From extended fig. 3, I can only see lower rate from Ks on single copy\n57 genes, and the rest look all the same. But in the text, it was stated that “B. insignis and\n58 its close relatives (Woodwardia prolifera and Pronephrium simplex) showed\n59 significantly lower dS and dN values than other ferns (Extended Data Fig. 3).” I just\n60 don’t see the pattern and what stats were done to support this?\n62 analysis. We agree that our original wording was too strong. The statement that “B.\n63 insignis has the lowest substitution rate among core leptosporangiate ferns” was an\n64 overgeneralization and did not adequately convey the nuances of the data. In reality,\n65 different genes evolve under different selective pressures, so it is unlikely that B.\n66 insignis would exhibit uniformly lower substitution rates across all genes compared to\n67 other ferns. Furthermore, our current sampling of fern species remains limited and may\n68 not fully capture the range of substitution rate variation in this clade.\n70 somewhat slower than many other core leptosporangiate fern: 1) WGD peak\n71 comparison: When analyzing the synonymous substitution (Ks) peaks corresponding to\n72 the shared ancient WGD event in ferns, B. insignis shows a significantly lower peak\n73 value (~1.7) compared to the fern Adiantum capillus-veneris (~2.25). This indicates\n74 slower accumulation of substitutions in B. insignis since that duplication event. 2)\n76 indicate slower molecular evolution in B. insignis in many comparisons. However, the\n77 differences are not statistically significant when B. insignis is compared to certain\n78 species with similarly slow rate (e.g., A. spinulosa, P. simplex, and B. hekouensis). We\n79 have now performed additional statistical tests and included them in Extended Data Fig.\n80 3 to quantify these differences and their significance. 3) Single-copy gene Ks values:\n82 to most other ferns examined. We consider this particularly informative because\n83 synonymous substitution rates (dS) are less subjected to selection and other\n84 confounding factors, providing more direct indicator of the neutral mutation rate. The\n85 consistently lower dS for B. insignis single-copy genes suggests a broadly slower base\n86 substitution rate.\n88 (Manuscript Lines 186-189). Additionally, we added a more comprehensive statistical\n89 comparison of dN and dS across species in Extended Data Fig. 3. This update provides\n90 a clearer visualization and statistical support to back the statement that B. insignis tends\n91 to have a slower evolutionary rate than many core leptosporangiate ferns (while\n92 acknowledging it is not the slowest in every comparison).\n94 expansions/contractions did not necessarily happen to B. insignis specifically, but the\n95 lineage leading to it. The closest species on this CAFÉ analysis was >100 million years\n96 apart. Therefore, there is no evidence that the expanded lignin biosynthesis genes “plays\n97 a crucial role in the persistence of this monotypic genus”\n99 contraction events identified by CAFÉ likely occurred in an ancestral lineage (over 100\n101 original wording in the manuscript attributing these changes specifically to B. insignis\n102 (and suggesting they “play a crucial role in the persistence of this monotypic genus”)\n103 was potentially misleading. We have revised that statement in the manuscript to be more\n104 precise and cautious, removing any implication that we have direct evidence for such a\n105 cause-and-effect in B. insignis.\n107 analysis with denser taxon sampling to shorten the phylogenetic branch lengths\n108 between B. insignis and the comparison species. By including additional closer relatives,\n109 we reduced the divergence times and thus minimized attributing the gene family\n110 changes to deep ancestral nodes. The updated results (included in Response Fig. 1 as\n111 referenced above) still indicate an expansion of lignin biosynthesis-related gene\n112 families on the branch leading to B. insignis. This gives us more confidence that the\n113 trend is relevant to B. insignis’s lineage. However, we now discuss these results with\n114 appropriate context: the expansion likely occurred in an ancestor and may have been\n115 retained in B. insignis due to its woody growth form. In the revised Manuscript, we\n116 have tempered our language about the significance of this finding for B. insignis’s\n117 persistence, noting it as a hypothesis rather than a confirmed mechanism.\n119 sounds like you have evidence for local adaptations. But I’m not sure if that evidence\n120 is presented.\n122 correct that our original phrasing was premature. Simply observing long-term isolation\n123 and genetic structure among the YN, SC, and VN lineages (from our population\n124 demographic analyses) does not, on its own, prove local adaptation. We agree that\n125 evidence for local adaptation should come from explicit tests, such as detecting\n126 selective sweeps or genotype–environment correlations, rather than from population\n127 structure alone.\n129 about “lineage-specific adaptations” that were not directly supported by evidence. More\n130 importantly, we strengthened our analysis to present comprehensive evidence\n131 supporting local adaptation through: 1) Selective sweep analysis: We conducted\n132 selective scans within each lineage (especially the SC and YN groups) and identified\n133 sets of genes that show signatures of recent positive selection unique to each lineage.\n135 biological sense in the context of their different habitats. This suggests that each lineage\n136 has experienced adaptive divergence; and 2) Genotype-environment association (GEA):\n138 = 0.4301, p < 0.05) and isolation-by-environment (Mantel’s r = 0.4058, p < 0.05)\n139 patterns in core adaptive loci. Redundancy analysis further demonstrated three distinct\n140 genetic clusters corresponding to each lineage (YN, SC, VN), with admixed\n141 populations occupying intermediate positions. Notably, the YN lineage’s genetic\n142 composition is strongly correlated with high diurnal temperature range (bio2), while\n144 precipitation in the driest month (bio14). The VN lineage formed a unique\n145 environmental cluster, with admixed populations showing heterogeneous\n146 environmental signatures.\n148 adaptation, with each lineage’s genome shaped by the specific environmental\n149 conditions of its region alongside the effects of geographic isolation. We have clarified\n150 the language to differentiate clearly between observed evidence of local adaptation (as\n151 above) and any remaining speculative inference. We believe this addresses the concern\n152 by replacing unsupported implication with concrete data.\n154 adaptability”.\n156 between “slow substitution rates” and “ecological adaptability” is rooted in the idea that\n157 some long-lived or “living fossil” plant species occupy exceptionally stable niches. For\n158 such species, a slow rate of molecular evolution is often a consequence of their enduring\n159 stability and success in an unchanging environment, rather than a cause of adaptability\n160 per se. In other words, if a species like B. insignis has persisted for millions of years in\n161 a relatively stable ecological niche, it may not require rapid genetic change to remain\n162 well-adapted. Thus, its substitution rate can appear low because it has already achieved\n163 a highly adapted state and continues to thrive without drastic genomic innovation. For\n164 example, Xiang et al. (2024) reported that the living fossil tree Euptelea pleiosperma\n165 exhibits an unusually low substitution rate, which the authors associated that specie’s\n166 long-term occupation of a conserved niche. Therefore, we suggest that in B. insignis,\n167 the observed slow substitution rate might be a signal that it is already well-adapted to\n168 its ecological niches and does not need rapid genetic turnover to cope with its\n169 environment. To avoid misunderstanding, we have rephrased this discussion in the\n170 manuscript. We now clarify that a slow substitution rate is interpreted as a sign of long-\n171 term persistence in stabile niches (Manuscript Line 391), and we avoid using the term\n172 “ecological adaptability” which was confusing.\n174 unusually slow rate of molecular evolution in the Cretaceous Eupteleaceae. Science\n177 about.\n179 a sudden spike in gene density on Chr04. We speculate that this localized increase in\n180 gene number is due to gene duplication events (possibly tandem duplications) or other\n181 adaptive processes that lead to a concentrated proliferation of related genes in this\n182 region.\n184 region underlying this spike. Our analysis localized the spike to a 1.2 Mb interval\n185 (Chr04:106.2-107.4 Mb) encompassing 91 annotated genes. Strikingly, many of these\n186 genes are closely related: our whole-genome collinearity analysis identified 79 gene\n187 pairs within this interval that reside in local syntenic blocks (indicative of recent\n188 duplications). Additionally, 32 of these genes in this cluster encode proteins with\n189 multiple conserved domains (e.g., Pfam:Photo_RC, Pfam:Oxidored_nitro), with many\n190 sharing identical or highly similar functional domains (Response Fig. 2). These findings\n191 strongly suggest that this sudden spike region on Chr04 results from tandem gene\n192 duplication events that expanded a family of gene with related functions.\n193\n196 assess if the sampling is representative and capturing the geography.\n198 distributed across tropical Asia, with documented occurrences spanning from 98.5°E to\n199 126.4°E in longitude and from 7.8°S to 25.5°N in latitude, according to GBIF (Global\n201 encompasses regions across Southeast Asia, extending from continental Asia (including\n202 known occurrences in China’s Yunnan) through the Indonesian archipelago. In this\n203 study, we sampled 94 individuals from 32 wild populations covering a broad portion of\n205 across Southeast Asia, with Vietnam’s representative population sampled at 12.19°N.\n207 insignis’ distribution. However, we would like to highlight the significance of this effort\n208 within the field of fern genomics. Given that B. insignis possesses an exceptionally\n209 large genome (~8 Gb)—ranking among the largest fern genomes ever sequenced at a\n210 population scale—this study represents one of the most comprehensive whole-genome\n211 resequencing initiatives conducted in ferns to date. We are confident that our sampling\n212 strategy captured the major geographic and environmental diversity of B. insignis\n213 across its core range, and thus provides robust insights into the species’ overall genetic\n214 structure and demographic history.\n216 known occurrences of B. insignis has been shown as Supplementary Fig. 9. As shown\n217 below:\n218\n220 insignis.\n222 records from NSII (National Specimen Information Infrastructure, China) and CVH\n223 (Chinese Virtual Herbarium), and GBIF (Global Biodiversity Information Facility)\n224 records, respectively.\n228 assembly and annotation of the Brainea insignis genome, as well as a population\n229 genomic assessment of the genetic diversity and demographic history of this\n230 endangered fern species. After several readings of the manuscript and supplemental\n231 information, I cannot find any major flaws in their study design, methods, or data\n232 analyses related to the genome sequencing, assembly, annotation, or populations\n233 genomic analyses. This work (with one exception) meets the standards of the field, and\n235 assessing their evaluation of genetic offset under future climate scenarios–this is too\n236 far afield from my expertise. I do, however, take some issue with some of their wording\n237 and interpretation of results, and I will detail a few examples at the end of this review.\n239 fern evolution and population genetics, and conservation genomics. In my opinion, the\n240 greatest impact of this publication is that it provides the first chromosome -level\n241 assembly for the Eupolypod II ferns, representing approximately 20% of extant fern\n242 diversity. It fills a major gap in genomic resources for ferns. However, Brainea insignis,\n243 a monotypic genus and “living fossil” akin to Gingko biloba and Cercidiphyllum\n244 japonicum, is an odd, atypical choice to represent this group. Perhaps, this choice\n245 reflects the authors’ focus on conservation, rather than comparative genomics for this\n246 study? Nonetheless, this genome assembly is a novel and needed resource for\n247 meaningful comparative genomics investigations of chromosome number, genome size,\n248 genic content, and rates of evolution relative to other ferns and to other major plant\n249 lineages, especially their sister group the seed plants. The authors briefly touch on this\n250 importance throughout the manuscript by confirming an ancestral WGD for the core\n251 leptosporangiate ferns, comparing intron sizes/exon numbers across ferns, illustrating\n252 the positive relationship between TE content and genome size in ferns, assessing rates\n253 of evolution among ferns, etc. In my opinion, the appeal and application of this\n254 manuscript would be broadened by refocusing and expanding discussion of these topics.\n256 the insightful suggestions for improving the manuscript. We also acknowledge that B.\n257 insignis is a monotypic genus and something of a “living fossil”, which indeed makes\n258 it an unusual representative of Eupolypod II. We chose B. insignis partly because of its\n259 conservation importance - it is endangered and thus of high conservation priority - and\n260 partly because of its pivotal phylogenetic position as one of the most divergent core\n261 leptosporangiate ferns. This dual rationale (conservation and comparative significance)\n262 is now clarified in the Introduction of the revised manuscript (Manuscript Lines 86-89).\n264 genomic stasis and divergence in ancient lineages. For example, our results support that\n265 the core leptosporangiate ferns (which include Eupolypod I and II clades, together ~60%\n266 of extant ferns) did not undergo any additional whole-genome duplications (WGDs)\n267 beyond the ancient WGD shared by all core leptosporangiates. Instead, other factors –\n268 such as the timing of lineage divergence and especially the accumulation of repetitive\n269 elements – appear to have had a greater influence on genome size variation among these\n270 ferns. We briefly touched on these points in the original manuscript, but we agree with\n271 the reviewer that they warrant a more prominent discussion.\n273 comparative analyses was somewhat limited (the closest relative genomes available,\n275 million years ago). We have proceeded with caution to avoid overinterpreting\n276 differences that could simply be due to deep evolutionary divergence. Nonetheless, the\n277 increasing availability of fern genomes will enable large-scale and in-depth\n278 comparative genomic analyses in the future, and such broader comparisons—including\n279 genome size evolution, TE dynamics, evolutionary rate variation, and synteny across\n280 ferns—would strengthen our study.\n282 structure in a biogeographic context, demographic history, gene flow, and signatures of\n283 local selection. The study of population genomics in ferns, particularly employing\n284 whole genome resequencing, is essentially non-existent and has been long hampered\n285 by a lack of genomic resources. I can think of only one other study (Yi et al. 2024. The\n287 estimate signatures of inbreeding and mutational load in a fern. I believe this manuscript\n288 will become a crucial resource/protocol for other researchers interested in population-\n289 level inquiries in ferns, in particular, for address many long-standing but previously\n290 intractable questions related to breeding systems in organisms that are spore dispersed\n291 and have independent sporophyte/gametophyte life phases.\n293 of B. insignis as a valuable resource for the community. As you pointed out, whole-\n294 genome resequencing studies in ferns are extremely scarce, primarily due to the\n295 historical lack of reference genomes and the computational challenges posed by the\n296 typically large fern genomes. In fact, B. insignis’s genome (~8 Gb) is among the largest\n297 fern genomes tackled with population resequencing to date, and this required us to\n298 overcome substantial technical and computational hurdles (e.g. handling massive\n299 amounts of sequencing data and repetitive sequences).\n301 to inspire and enable more studies on fern population genetics—particularly to address\n302 the unique questions you mentioned, like how breeding systems and life-cycle\n303 characteristics of ferns influence their genetic diversity and inbreeding levels.\n305 evolution of lignin syringyl “S” units, mainly reported from angiosperms, in B. insignis.\n307 2024. https://doi.org/10.1038/s41477-025-684 01978-y), indicating that lignin S-units\n308 have evolved repeatedly across land plants, and that ferns have independently evolved\n309 mechanisms of structural support and vasculature in terrestrial environments.\n311 our results suggest that B. insignis (and possibly other ferns) have unique expansions\n312 in certain lignin-biosynthesis gene families, which could relate to the production of\n313 lignin S-units or analogous compounds for mechanical support. The recent finding by\n315 units—once thought to be exclusive to angiosperms—have evolved independently in\n316 multiple plant lineages. This independent evolution of lignin chemistry and vascular\n317 support strategies in ferns versus seed plants provides important context for our\n318 observations. The evolution of lignin S-units in ferns remains an intriguing open\n319 question. Our study not only supports the idea of independent lignin pathway evolution\n320 but also sets the stage for future research to explore how ferns achieve their structural\n321 support.\n323 inaccurately depicted background information or have interpreted their results too\n324 liberally. Below are three examples that could be amended:\n326 bearing vascular plants…” This may seem overly picky, but I think it misrepresents the\n327 placement/informativeness of ferns in plant evolution. Ferns do not “bridge” or have\n328 some type of “intermediate” phylogenetic placement. (As the authors know) ferns are\n329 the sister group to seed plants, and share a common ancestor with seed plants, not shared\n330 with bryophytes. Please clarify for accuracy.\n332 our original phrasing was inaccurate. Ferns are nor an intermediate between bryophytes\n333 and seed plants; rather, ferns are the sister group to seed plants. This means ferns and\n334 seed plants share a common ancestor that is not shared with bryophytes. We have\n335 corrected this sentence in the Introduction of the revised manuscript to accurately\n336 reflect fern phylogeny. The revised statement is: “As the sister group to seed plants,\n337 ferns share a common ancestor with them, and thus provide crucial insights into plant\n338 evolution, particularly transitions in anatomy, life history, and reproductive strategies”.\n340 2-Line 130 “These results collectively underscore the high contiguity, consistency, and\n341 completeness of our B. insignis assembly.” I am concerned about the authors’\n342 conclusion of completeness given the BUSCO score of 81.5% complete recovery (77.7%\n343 based on annotation assessment, Supplementary Table 3). Despite having high LAI and\n345 usually 95% complete or higher). At minimum, the authors need to address the possible\n346 reasons for the low BUSCO recovery and discuss whether it is an artifact of their\n347 assembly or possibly reflecting a biological reality (unlikely).\n349 regarding assembly “completeness.” In the original manuscript, our BUSCO analysis\n350 reported approximately 81.5% completeness for the genome assembly (and about 77.7%\n351 for the annotated gene set), which is indeed lower than the typical >90–95% seen in\n352 many seed plant genomes. We agree that we needed to provide an explanation and\n353 adjust our claims accordingly.\n355 along with the lack of a high-quality fern reference genome. As mentioned in our\n356 response to Reviewer #1, the use of the embryophyta_odb10 dataset, which is biased\n357 toward angiosperms, may not have been optimal for our study. In the revised version,\n358 we repeated the BUSCO analysis with the updated viridiplantae_odb12 dataset\n359 (updated July 1, 2025), which better represents non-seed plants. This more\n360 comprehensive BUSCO set greatly improved our completeness estimates. The genome\n361 assembly now scores 97.4% Complete, and the gene annotation scores 84.2%, as shown\n362 in the updated Supplementary Table 1 (and Response Table 1). The large jump in the\n363 assembly’s BUSCO score (from ~81% to ~97%) indicates that the prior lower score\n364 was largely an artifact of using an angiosperm-skewed BUSCO dataset\n365 (embryophyta_odb10) on a fern genome. In other words, many of the “missing”\n367 divergence or absence from the embryophyta BUSCO list.\n369 viridiplantae_odb12 dataset\nEvaluation Complete Single-copy Duplicated Fragmented Missing Total\nType (C) (S) (D) (F) (M) BUSCOs (n)\nGenome\n97.4% 89.9% 7.5% 2.1% 0.5% 822\nAssembly\nGene\n84.2% 79.1% 5.1% 12.7% 3.2% 822\nAnnotation\n370 3-Line 187 “relative rate tests suggest that B. insignis has the slowest evolutionary rate\n371 among core leptosporangiate ferns (Supplementary Table 7).” Given the information in\n373 have a significantly slower rate of evolution (in 31 single copy genes) relative to several\n374 other core leptosporangiate ferns, including L. chinensis and D. repens. Please amend\n375 this statement and temper your conclusion that “this slow evolutionary rate may have\n376 bolstered genomic stability, enabling…”.\n378 statements on substitution rates. This comment echoes the concern raised by Reviewer\n379 #1, and we have taken it very seriously. We have amended the manuscript to remove\n380 the claim that B. insignis has the single “slowest” rate among core leptosporangiate\n381 ferns and to significantly soften the interpretation about what a slower rate implies for\n382 the species.\n384 recognize that our original wording was too strong and could be misleading. We have\n385 rewritten the relevant passage in the Results/Discussion to clarify that while B. insignis\n386 shows generally slower substitution rates compared to many ferns, it is not an outlier in\n387 all comparisons. For instance, as you pointed out, species like Lygodium chinensis and\n389 insignis in our single-copy gene analyses."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I believe the authors have nicely addressed most of my previous concerns. I found the new synteny comparison with\nAlsophila very interesting. I do want to point out that computing “Ks” between LTR repeats is not valid. Ks is synonymous\nsubstitution and would only apply to protein-coding genes. How can you calculate Ks for repeats?", "reviewer_2": "Dear Editors and Authors,\nIt is my opinion that Li et al. have satisfactorily addressed both reviewers’ comments in a manner that improves the execution\nof the study and better highlights its immediate and broader impacts. The methodologies are sound, of the standards of the\nfield, reproducible, and noteworthy. As described in my earlier review, I believe this manuscript is a well-written and valuable\ncontribution as the first complete Eupolyplod II fern genome assembly and as a rare example of using whole genome\nsequencing to investigate population level diversity for a fern species. I especially appreciate that the authors have\nbroadened their discussion of the importance of this genome for investigating fern genome evolution. Of particular note, is\ntheir result that large genome size in the polypods seems largely dictated by TE expansion, rather than by WGD followed by\nslow rates of diploidization – this is important because it counters the existing paradigm for how ferns have accrued large\nchromosome numbers and genome sizes relative to angiosperms. I believe this study will also serve as a useful and well-\ncited model for population genomics approaches in ferns, a group historically underrepresented in population genomics\nbecause of limited genomic resources and large population size.\nOne caveat is that I wish the manuscript had been assigned to a reviewer with expertise in the evaluation of genetic offset\nunder future climate scenarios. I feel that this aspect of the manuscript largely went unevaluated and might benefit from an\nexpert’s opinion.", "reviewer_3": "Xia et al. present a study that assembled a de novo genome of a cycad fern and conducted population genomic analyses\nbased on 94 individuals. Although the dataset is substantial, primarily due to the species’ large genome size, I find the\nnovelty of the analysis somewhat underwhelming. Below, I outline several concerns and suggestions that could enhance the\nstudy’s relevance and originality, particularly in comparison to other population genomic studies on species with smaller\ngenome sizes.\nMain Comments\n1.Importance of the Genome and Comparative Analysis: While the de novo genome assembly of this endangered species\nwithin the Eupolypods II clade is significant, the genomic analyses presented in the study feel overly general. Beyond the\nbasic genomic structure in Figure 1, Figure 2 highlights WGD (whole-genome duplication) events and phylogenetic\nrelationships across diverse species. It would be more impactful to compare this genome with other published fern genomes\nto examine genome collinearity and divergence among closely related species. A more detailed analysis could investigate\nunique transposable elements (TEs) or gene expansions/contractions specific to this species, as well as large-scale\nstructural variations in the genome. Linking these variations to associated genes and their functional enrichment would\nprovide deeper insight into this species’ unique genomic features.\n2.Integration of Demographic, Divergence, and Gene Flow Models: The study constructs demographic, differentiation, and\ngene flow models separately. However, these processes are interrelated and should ideally be integrated into a cohesive\nframework. Typically, divergence models incorporate demographic and gene flow histories to provide more accurate and\ncomprehensive inferences. Analyzing these processes separately can lead to confounding effects and compromise the\naccuracy of the results. I suggest the authors revisit their fastsimcoal analysis, integrating demographic and gene flow\nscenarios into unified models for better reliability and interpretability.\n3.Genetic Load and Population Decline: One of the study's key contributions is its emphasis on the endangered status of this\nspecies and the impact of population decline on genetic variation and genetic load. However, the current analysis only\ncompares genetic load across the three genetic groups, which exhibit minimal differences and provide limited insight into the\nextent of genetic load accumulation. To enhance this analysis, the authors could employ established methods for measuring\ngenetic load, such as the 0-fold/4-fold ratio or the homozygous/heterozygous counts of loss-of-functon(LoF) and deleterious\nmutations. Furthermore, comparing the genetic load of this species to data from other studies on species with similarly small\npopulation sizes or those in an endangered state would offer a more meaningful context. This comparison could help clarify\nthe severity of the genetic load accumulation in this fern species and its relative position on the spectrum of endangerment.\n4.Issues with the Ratio of Nonsynonymous to Synonymous Diversity: On line 312, the authors report a nonsynonymous-to-\nsynonymous diversity ratio close to 1, which suggests a potential calculation error. I recommend carefully reviewing the\nmethodology used to compute this ratio. If the result is accurate, it represents an unusual or unique finding, as such a high\nvalue (π0-fold/π4-fold) is rarely observed. In that case, further investigation is necessary to determine the underlying causes\nof this pattern.\n5.Genetic Vulnerability to Climate Change: The study includes an analysis of genetic vulnerability to future climate change,\nbut this is now a fairly standard approach. It’s important to critically assess whether such methods are suitable for this\nspecies and whether they provide accurate predictions. The genetic offset analysis assumes environmental adaptation, but\nin this case, geographic distance is strongly correlated with environmental distance. This makes it difficult to disentangle the\neffects of selection-driven environmental adaptation from neutral genetic drift due to geographic separation. To address this\nissue, the authors could focus on specific identified adaptive variants and test whether they are associated with\nenvironmental variables, geographic differences, or a combination of both. Differentiating these factors is essential for\nunderstanding the role of environmental adaptation before conducting genetic offset analyses.\n6.Quantification of Genetic Offset Results: The genetic offset analysis currently presents results in the form of a colored plot,\nwhich makes it difficult to pinpoint which populations are more vulnerable. I recommend quantifying the offset values (e.g.,\nlocal, forward, and reverse offsets) for each population or genetic cluster and visualizing them using boxplots or similar\napproaches. This would allow readers to clearly identify which populations are better adapted and which are more\nvulnerable to future climate change.\nMinor Comments\n1.On line 241, the term core variants is unclear. How is it defined? Providing a clear definition would improve clarity.\n2.In general, the Results section contains several ambiguous statements, partly because the methods are described later.\nRestructuring or clarifying these sections would make the study easier to follow."}, "author_response": "Responses to Reviewers’ comments\nReviewer #1 (Remarks to the Author):\nLine 184-186. First, I won’t consider Pronephrium simplex a “close relative”. It’s in a\ndifferent family. Second, I still don’t understand why only single copy genes show rate\ndifference, and not the low copy genes. Importantly, since single copy genes only\naccount for a small proportion of the genome, can they really represent the evolutionary\nrate of whole genomes?\nResponse: We thank the reviewer for the taxonomic clarification and agree that\nPronephrium simplex is not a close relative of Brainea insignis, we have therefore\nremoved this phrasing from the manuscript.\nWith respect to the use of only single-copy genes, we acknowledge that they represent\nonly a small fraction of the genome and thus may not fully capture genome-wide\nsubstitution rates. Nevertheless, analysis of low-copy or multi-copy genes are more\nvulnerable to confounding by complex gene family dynamics, including duplication\nand loss, sub- or neofunctionalization, and episodes of positive selection in individual\nparalogs. These processes introduce rate heterogeneity that can obscure lineage-level\nsignals. For this reason, strict single-copy orthologs remain the standard approach for\nphylogenomic and molecular rate comparisons, as they minimize the risk of paralogs\nand other biases.\nWe also recognize the limitations of this approach, particularly given the large\nevolutionary distances among available fern genomes, and the fact that most taxa are\nrepresented by transcriptomes rather than high-quality whole genomes. To provide an\nindependent line of evidence, we therefore conducted a genome-wide synteny analysis\nwith the tree fern Alsophila spinulosa, a lineage previously reported to evolve slowly.\nDespite divergence over 200 Mya (TimeTree5), we detected extensive 1:1 syntenic\nblocks between B. insignis and A. spinulosa (Supplementary Fig. 6). This level of\nconserved collinearity supports our interpretation that B. insignis genome has evolved\nat a relatively slow pace, at least in terms of genomic architecture. This result has been\nadded in the revised manuscript (Lines 189-193) and is presented in the supplementary\nmaterial.\nTogether, these complementary analyses strengthen the inference of slow genome\nevolution in B. insignis while underscoring the need for additional high-quality fern\ngenomes to enable broader and more precise assessments of evolutionary rate variation.\nSupplementary Fig. 6 | Synteny analysis between Alsophila spinulosa and Brainea\ninsignis.\nLine 189-192. The relative rate tests similarly have issues. I am not sure to what extent\n31 single copy genes are representative of this species. Further, with Psilotum nudum\nas the outgroup, the timescale here is vast. In other words, is the rate difference truly\nspecific to Brainea or the hundreds of million years of evolution leading to it? I am not\ncertain the analyses here are valid.\nResponse: We appreciate the reviewer’s valuable comments. We agree that a set of 31\nsingle-copy genes is limited and, on its own, cannot decisively represent genome-wide\nsubstitution rates. Accordingly, we have supplemented our study with results from\ngenome synteny analyses, which clearly show that, compared with the more distantly\ndiverged Alsophila, Brainea has retained a substantial amount of syntenic blocks\n(Supplementary Fig. 6). These findings provide structural evidence that Brainea has\nindeed evolved at a relatively slower pace, at least in terms of preserving genomic\narchitecture.\nRegarding the relative-rate test, we applied the classical pairwise approach of Tajima,\nwhich compares the substitution counts between ingroup taxa relative to an outgroup O\nto test whether A–O and B–O distances differ significantly. We agree that relying solely\non a distantly related outgroup such as Psilotum nudum may introduce artifacts due to\nsubstitution saturation over deep timescales. To mitigate this, we repeated all tests with\ntwo additional outgroups spanning shallower and intermediate divergence—Marsilea\nvestita and Asplenium formosae. As summarized in Response Table 1, the majority of\nsignificant comparisons consistently indicate a slower rate in B. insignis across\noutgroups. Concordance of the signal across multiple outgroup depths increases\nconfidence that the effect reflects the Brainea lineage rather than deep stem\naccumulation.\nTaken together—(i) the outgroup-robust relative-rate results and (ii) the independent\nevidence of conserved large-scale collinearity—our inference is that B. insignis has\nevolved comparatively slowly, at least in terms of genomic architecture, while we\nexplicitly acknowledge the limits imposed by current fern genomic resources and the\nsmall single-copy gene set. Further high-quality genomes from closely related lineages\nwill enable more precise, genome-wide rate estimation.\nResponse Table 1 | Relative rate tests comparing B. insignis with other ferns using\nM. vestita and A. formosae as outgroups separately\nIngroup1 Ingroup2 Outgroup Genes Identical Divergent Ingroup1 specific Ingroup2 specific Outgroup specific Slow CHI^2_test_statistic P-value\nAlsophila spinulosa B. insignis M. vestita 31 9261 694 567 628 1310 - 3.11 0.07763\nAdiantum capillus-venerisB. insignis M. vestita 31 9426 714 602 457 1499 B. insignis 19.85 0.00001\nAsplenium formosae B. insignis M. vestita 31 9273 558 488 351 1719 B. insignis 22.37 0.000001\nWoodwardia prolifera B. insignis M. vestita 31 10017 231 143 189 2361 W. prolifera 6.37 0.01158\nPronephrium simplex B. insignis M. vestita 31 9773 313 241 246 2174 - 0.72 0.588\nBolbitis hekouensis B. insignis M. vestita 31 8905 429 304 271 1832 B. insignis 1.89 0.16876\nDavallia repens B. insignis M. vestita 31 8993 496 377 314 1754 B. insignis 5.74 0.01655\nLoxogramme chinensis B. insignis M. vestita 31 7990 511 441 309 1473 B. insignis 23.23 0.000001\nNephrolepis cordifolia B. insignis M. vestita 31 9282 429 268 302 1845 - 2.03 0.15442\nIngroup1 Ingroup2 Outgroup Genes Identical Divergent Ingroup1 specific Ingroup2 specific Outgroup specific Slow CHI^2_test_statistic P-value\nWoodwardia prolifera B. insignisA. formosae 31 10968 124 197 242 1079 W. prolifera 4.61 0.03173\nPronephrium simplex B. insignisA. formosae 31 10703 166 325 338 917 - 1.31 0.28\nBolbitis hekouensis B. insignisA. formosae 31 9987 243 436 357 790 B. insignis 7.87 0.00503\nDavallia repens B. insignisA. formosae 31 10125 273 565 376 759 B. insignis 37.96 0.000001\nLoxogramme chinensis B. insignisA. formosae 31 8929 274 662 354 658 B. insignis 93.37 0.000001\nNephrolepis cordifolia B. insignisA. formosae 31 10454 234 422 364 823 B. insignis 4.28 0.03857\nLine 200-201. “Larger fern genomes tend to harbor a higher proportion of LTR-RTs\nwith earlier insertion time”. How can you accurately estimate and compare insertion\ntimes when there is rate heterogeneity among fern lineages (a point you are trying to\nmake above)?\nResponse: We agree that lineage-specific rate heterogeneity complicates cross-species\ncomparisons of LTR-RT insertion ages, which are typically inferred from divergence\nbetween the 5′ and 3′ LTRs.\nTo address this variation, we adopted a phylogenetically informed, relative calibration\nstrategy. Specifically, we first estimated synonymous substitution rates (Ks) from\northologous genes across fern lineages, using fossil-calibrated divergence times from\nthe TimeTree5 database to derive lineage-specific substitution rates (r = Ks_ortholog /\n2t). These rates were then used to rescale observed Ks values, making them comparable\nwith the paranome Ks distribution of the focal species (via wgd-v2, developed by our\ncollaborator Yves Van de Peer’s group). This normalization corrects for background\nmutation rate differences, enabling cross-species comparisons of Ks distributions and\nalignment with shared evolutionary events such as WGDs. For LTR retrotransposons,\nwe used LTR_retriever to calculate Ks between the 5’ and 3’ LTRs of each element.\nRather than interpreting raw LTR-Ks values directly, we focused on relative patterns\nwithin the rescaled ortholog-based framework, which allowed us to assess shifts in the\nmode or peak of LTR-Ks distributions in a phylogenetically consistent manner. The\nobservation that species with larger genomes retain not only more LTR-RTs but also\nproportionally older insertions—even after correction—supports our conclusion that\ngenome size expansion is linked to the persistence of ancient LTR-RTs. While we\nrecognize this represents a first-order correction and some residual rate variation may\nremain, the consistency of the results across lineages gives us confidence in the\nrobustness and biological significance of this signal.\nIn summary, while absolute insertion times should be interpreted cautiously, our use of\ntime-calibrated orthologous divergence to normalize rates provides a phylogenetically\njustified framework for comparing the relative ages of LTR-RT bursts across ferns. To\nenhance transparency and reproducibility, we added a detailed description of this rate\nnormalization procedure in the Methods section (Lines 494-498).\nLine 213-216. The vast majority of the genes are evolving under purifying selection\n(w<1) and I do not believe this can be used to support “strong functional constraint on\nthese lignin-related genes across these fern lineages.” In addition, Brainea and\nSphaeropteris span deep evolutionary time, during which arborescence only accounts\nfor a small blip in time – does this calculation actually mean anything?\nResponse: We thank the reviewer for raising this important point. We agree that ω < 1\nis a general pattern for most protein-coding genes, and we also acknowledge that the\ndivergence between Brainea and Sphaeropteris spans a deep evolutionary timescale,\nduring which the arborescent growth form represents only a relatively small fraction.\nTo further clarify, we extracted the lignin-related genes identified in B. insignis and\nexamined their micro-synteny blocks with Sphaeropteris lepifera and Alsophila\nspinulosa. The results (presented in Supplementary Fig. 10) show that a subset of\nlignin-related genes is highly conserved across these species. Thus, our intention is not\nto imply a direct link between ω ratios and arborescence, but rather to highlight the\nconserved nature of lignin pathway genes, consistent with their central role in cell-wall\nbiology across ferns. We have therefore softened our original interpretation and no\nlonger present ω values as primary evidence. Instead, we highlight an independent line\nof support: conserved microsynteny among a subset of lignin-pathway genes in Brainea\ninsignis, Sphaeropteris lepifera, and Alsophila spinulosa (Supplementary Fig. 10). The\nlong-term preservation of gene order and content across these deeply divergent lineages\nis more consistent with functional and dosage constraints acting on pathway\narchitecture than with ω values alone, which can be ambiguous or misleading.\nSupplementary Fig. 10 | Microsynteny analysis of lignin-related genes across\nthree fern species\nLine 219-222. Another possibility is that the expansion of these gene families has\nnothing to do with woodiness!\nResponse：We agree that CAFÉ-inferred expansions do not, by themselves, establish\nany causal relationship with woodiness. CAFE detects net turnover in family size along\nbranches under a birth–death model; it is agnostic to function and trait causation. To\navoid over-interpretation, we have removed the speculative causal statement and now\npresent these results as descriptive patterns of gene-family turnover. We explicitly note\nthat linking expansions to woodiness would require additional evidence—e.g.,\nreplicated associations across independent transitions using phylogenetic comparative\nmodels that control for genome size/WGD and life history, together with tissue-specific\nexpression and functional assays. The manuscript has been revised accordingly.\nLine 225-228. What’s the evidence for “convergent evolution of lignin specialization”?\nDo we know what types of lignin does Brainea actually have? There are a lot of vague\nassociations and assumptions in this section.\nResponse：We thank the reviewer for this valuable comment. We agree that our original\nphrase “convergent evolution of lignin specialization” exceeded the available evidence.\nIndeed, the exact lignin composition of B. insignis has not been characterized, and thus\ndirect support for convergent evolution is lacking. Our initial association between\nBrainea and Stenochlaena palustris was motivated by their shared family affiliation,\nwhich led us to consider potential similarities in lignin-related traits. To avoid over-\ninterpretation, we have removed this speculative conclusion and revised the text to\npresent the interpretation more cautiously (Lines 224-229).\nMinor comments:\nLine 60. “Ferns share a common ancestor with them” – Every organism on this planet\nshares a common ancestor with another, so this statement is not needed.\nResponse: Thank you for pointing this out. We acknowledge that the statement could\nbe seen as trivial given the universal common ancestry of life. It has been removed to\nimprove precision.\nLine 89. Who is “them” in “the broader biodiversity that depends on them”?\nResponse：We appreciate the reviewer’s meticulous comment. We have revised this\nsentence to: Given the species’ dual importance to conservation and evolutionary\nstudies, understanding its genomic makeup and evolutionary resilience is critical for\nguiding the conservation of both the species itself and the biodiverse habitats it sustains.\nLine 127. “Clipping-based Revealing Assembly Quality” should be “Clipping Reveals\nAssembly Quality” or “Clipping information for Revealing Assembly Quality”?\nResponse: We thank the reviewer for this suggestion. The term has been revised to\n“Clipping Information for Revealing Assembly Quality (CRAQ)” for improved clarity\nand accuracy.\nReviewer #2 (Remarks to the Author):\nDear Editors and Authors,\nIt is my opinion that Li et al. have satisfactorily addressed both reviewers’ comments\nin a manner that improves the execution of the study and better highlights its immediate\nand broader impacts. The methodologies are sound, of the standards of the field,\nreproducible, and noteworthy. As described in my earlier review, I believe this\nmanuscript is a well-written and valuable contribution as the first complete Eupolyplod\nII fern genome assembly and as a rare example of using whole genome sequencing to\ninvestigate population level diversity for a fern species. I especially appreciate that the\nauthors have broadened their discussion of the importance of this genome for\ninvestigating fern genome evolution. Of particular note, is their result that large genome\nsize in the polypods seems largely dictated by TE expansion, rather than by WGD\nfollowed by slow rates of diploidization – this is important because it counters the\nexisting paradigm for how ferns have accrued large chromosome numbers and genome\nsizes relative to angiosperms. I believe this study will also serve as a useful and well-\ncited model for population genomics approaches in ferns, a group historically\nunderrepresented in population genomics because of limited genomic resources and\nlarge population size.\nOne caveat is that I wish the manuscript had been assigned to a reviewer with expertise\nin the evaluation of genetic offset under future climate scenarios. I feel that this aspect\nof the manuscript largely went unevaluated and might benefit from an expert’s opinion.\nResponse: We thank the reviewer for the thoughtful evaluation and endorsement. We\nhave clarified how the B. insignis genome informs fern genome evolution (e.g., WGD\ncontext, TE dynamics) and explicitly noted that, to our knowledge, this is the first\napplication of a genetic-offset framework in a fern. We also acknowledge the reviewer’s\ncaveat and, in the Discussion/Methods, more clearly justify the suitability and\nlimitations of genetic-offset analyses for Brainea (long-lived, perennial, spore-\ndispersed) and the steps we took to mitigate confounding (see also responses to\nReviewer #3, point 5).\nReviewer #3 (Remarks to the Author):\nXia et al. present a study that assembled a de novo genome of a cycad fern and\nconducted population genomic analyses based on 94 individuals. Although the dataset\nis substantial, primarily due to the species’ large genome size, I find the novelty of the\nanalysis somewhat underwhelming. Below, I outline several concerns and suggestions\nthat could enhance the study’s relevance and originality, particularly in comparison to\nother population genomic studies on species with smaller genome sizes.\nResponse: We appreciate the reviewer’s careful assessment and constructive\nsuggestions. Below we describe concrete changes and clarifications made to enhance\noriginality, interpretability, and rigor.\n1. Importance of the Genome and Comparative Analysis: While the de novo genome\nassembly of this endangered species within the Eupolypods II clade is significant, the\ngenomic analyses presented in the study feel overly general. Beyond the basic genomic\nstructure in Figure 1, Figure 2 highlights WGD (whole-genome duplication) events and\nphylogenetic relationships across diverse species. It would be more impactful to\ncompare this genome with other published fern genomes to examine genome\ncollinearity and divergence among closely related species. A more detailed analysis\ncould investigate unique transposable elements (TEs) or gene expansions/contractions\nspecific to this species, as well as large-scale structural variations in the genome.\nLinking these variations to associated genes and their functional enrichment would\nprovide deeper insight into this species’ unique genomic features.\nResponse: We appreciate the reviewer’s insightful suggestions. We agree that\ncomparative genomic analyses, such as collinearity assessments, investigation of\nspecies-specific transposable elements, gene family dynamics, and structural variation,\nprovide the most incisive route to uncover lineage-specific features. The current\nconstraint is data availability: chromosome-level fern genomes remain scarce, and\nwithin eupolypods II our B. insignis assembly is, to our knowledge, the first at this\nresolution. The nearest available references (Ceratopteris richardii, Adiantum capillus-\nveneris) diverged from Brainea over 100 million years ago and exhibit substantial\ndifferences in chromosome number, reflecting deep karyotypic reorganization. Coupled\nwith extensive gene turnover, large-scale rearrangements, and TE flux, this severely\nlimits the resolution of fine-scale “one-to-one” comparative analyses. To avoid over-\ninterpreting differences driven mainly by deep evolutionary divergence rather than true\nlineage-specific changes, we interpreted our comparative analyses with caution. We\nagree that as more high-quality fern genomes become available, broader and more\ndetailed comparative studies—including genome size evolution, TE dynamics, synteny,\nand gene family turnover across ferns—will be both possible and highly informative.\nIn the revised manuscript, we therefore performed genome-wide collinearity analysis\nbetween B. insignis and Alsophila spinulosa (diverged >200 Mya) and observed\nextensive 1:1 syntenic blocks (Supplementary Fig. 6), highlighting the structural\nconservation of the B. insignis genome. This analysis also revealed a shared WGD event,\nbut indicated that WGD alone cannot fully account for genome size variation in ferns;\nrather, LTR insertions and accumulation appear to play a more dominant role.\nFurthermore, micro-synteny analysis revealed that lignin biosynthesis–related gene sets\nare conserved between B. insignis, Sphaeropteris lepifera, and Alsophila spinulosa\n(Supplementary Fig. 10), despite their divergence over 200 Mya, underscoring the long-\nterm evolutionary stability of these functionally critical pathways. Taken together, these\nfindings provide important insights into fern genome evolution. We hope the reviewer\nappreciates the current limitations in fern genomic resources and recognizes that,\nalthough some analyses are standard in seed plants or animals, our study provides\npioneering insights into fern genomics and population genetics.\n2. Integration of Demographic, Divergence, and Gene Flow Models: The study\nconstructs demographic, differentiation, and gene flow models separately. However,\nthese processes are interrelated and should ideally be integrated into a cohesive\nframework. Typically, divergence models incorporate demographic and gene flow\nhistories to provide more accurate and comprehensive inferences. Analyzing these\nprocesses separately can lead to confounding effects and compromise the accuracy of\nthe results. I suggest the authors revisit their fastsimcoal analysis, integrating\ndemographic and gene flow scenarios into unified models for better reliability and\ninterpretability.\nResponse: We thank the reviewer for this important suggestion. We agree that\ndivergence, demography, and gene flow are inherently interdependent and ideally\nshould be inferred within a unified modeling framework. While fastsimcoal2 does\nsupport such joint inference, simultaneous estimation of many parameters from the site\nfrequency spectrum (SFS) can lead to challenges in parameter identifiability—such as\ntrade-offs between divergence time, effective population size changes, and migration\nrates—as well as flat likelihood surfaces, which may compromise the reliability of\nestimates. To balance realism with statistical robustness, we adopted a stepwise,\nhierarchical modeling strategy: we first fitted simpler models with fewer parameters\n(e.g., divergence without migration), then incrementally introduced additional\nprocesses (asymmetric continuous migration and pulse admixture), retaining\ncomplexity only when supported by the data. Practically, each step was evaluated by\nrepeated optimizations from diverse starting points, and model choice was based on\nAkaike Information Criterion (AIC) and ΔAIC. The best-supported models therefore\nintegrate divergence, demographic history, and gene flow rather than treating them in\nisolation, and their parameter estimates are reported accordingly (Supplementary Figs.\n17–19; Tables 8–10). Where appropriate, we also compared key time-scale parameters\nwith independent skyline-based inferences to assess concordance.\nSuch a tiered modeling framework has also been applied in recent studies—for example,\nXiao et al., (2023, PNAS) used a similar progressive approach to robustly reconstruct\nthe domestication history of grapevine. We have revised the text to make this workflow\nexplicit and to avoid the impression that processes were analyzed independently. The\nstepwise nesting serves only to ensure parameter identifiability; the final inferences are\ndrawn from the highest-supported integrated models (Lines 557–566).\nReference: Xiao et al., Adaptive and maladaptive introgression in grapevine\ndomestication. Proc. Natl. Acad. Sci. 24 (120), 2023.\n3. Genetic Load and Population Decline: One of the study’s key contributions is its\nemphasis on the endangered status of this species and the impact of population decline\non genetic variation and genetic load. However, the current analysis only compares\ngenetic load across the three genetic groups, which exhibit minimal differences and\nprovide limited insight into the extent of genetic load accumulation. To enhance this\nanalysis, the authors could employ established methods for measuring genetic load,\nsuch as the 0-fold/4-fold ratio or the homozygous/heterozygous counts of loss-of-\nfuncton(LoF) and deleterious mutations. Furthermore, comparing the genetic load of\nthis species to data from other studies on species with similarly small population sizes\nor those in an endangered state would offer a more meaningful context. This\ncomparison could help clarify the severity of the genetic load accumulation in this fern\nspecies and its relative position on the spectrum of endangerment.\nResponse: We agree that assessing genetic load with site-class–based metrics and\nhomozygosity of deleterious variants provide a clearer picture than groupwise\ncomparisons alone. Accordingly, we expanded the analysis in three ways. First, we\nderived 0-fold and 4-fold degenerate sites from the reference CDS and computed the\nratio of nonsynonymous to synonymous diversity (π₀/π₄) to gauge the efficacy of\npurifying selection (Supplementary Table 11). Second, we quantified the burden of\ndeleterious alleles by distinguishing homozygous and heterozygous genotypes for\npredicted loss-of-function (LoF) and deleterious variants, and we summarized the\nproportion of homozygous deleterious genotypes per individual as a proxy for realized\nload (Fig. 4e–f). Third, we interpreted these patterns in the context of recent population\ndecline, where increased drift and inbreeding are expected to elevate realized\n(homozygous) load even if the total number of deleterious alleles remains similar across\ngroups.\nFor cross-species context, resources in ferns are still limited; to our knowledge,\nAlsophila spinulosa is the only species with both a high-quality genome and population\nresequencing suitable for comparison. We therefore discuss Brainea insignis relative to\nA. spinulosa in the manuscript (Lines 327–337), noting life-history and reproductive\ndifferences that preclude direct quantitative benchmarking against seed plants. This\nprovides the most biologically relevant frame of reference currently available while\navoiding potentially misleading comparisons across very different clades.\nTogether, these additions move beyond coarse group comparisons to a site-class and\ngenotype-state perspective, offering a more informative view of genetic-load\naccumulation associated with recent demographic collapse in B. insignis, while clearly\nstating the comparative limits imposed by current fern genomic resources.\n4. Issues with the Ratio of Nonsynonymous to Synonymous Diversity: On line 312, the\nauthors report a nonsynonymous-to-synonymous diversity ratio close to 1, which\nsuggests a potential calculation error. I recommend carefully reviewing the\nmethodology used to compute this ratio. If the result is accurate, it represents an unusual\nor unique finding, as such a high value (π0-fold/π4-fold) is rarely observed. In that case,\nfurther investigation is necessary to determine the underlying causes of this pattern.\nResponse: We thank the reviewer for identifying this issue. The initially reported π₀/π₄\nvalue was inflated due to a mistake in the definition of degenerate sites: non-coding\nregions were inadvertently included in the annotation of 0-fold and 4-fold degenerate\nsites, which should be restricted to coding sequences (CDS). This contamination\ncompromised the accuracy of the ratio, rendering it unreliable. We have corrected the\nanalysis using properly filtered sites, and all results have been updated accordingly\n(Supplementary Table 11).\nWe re-ran the analysis using a codon-aware, population-aware pipeline consistent with\nrecent practice. Specifically, we (i) annotated true 0-fold and 4-fold sites from the\nreference CDS using a degenerate-site annotator (i.e., degenotate; Mirchandani et al.,\n2024), (ii) intersected these coordinates with the population VCF to retain only\npolymorphisms at those sites, and (iii) recomputed nucleotide diversity (π) with\nANGSD. The corrected π₀/π₄ estimates are substantially lower than the original near-\nunity value (see Supplementary Table 11), confirming that the initial result was an\nartifact of site misclassification.\nInterpreted biologically, the corrected ratios remain elevated relative to many published\nplant datasets (e.g., Plomion et al., 2018), which is consistent with relaxed purifying\nselection and increased drift/inbreeding associated with recent population decline. We\nhave revised the Results and Methods accordingly (Lines 314–316; 587–592) and are\ngrateful to the reviewer for prompting this correction.\nReference:\n[1] Mirchandani et al., 2024. High-throughput variant calling workflow for population\ngenomics. Molecular Biology and Evolution 41(1).\n[2] Plomion et al., 2018. Oak genome reveals facets of long lifespan. Nature Plants 4,\n440–452.\n5. Genetic Vulnerability to Climate Change: The study includes an analysis of genetic\nvulnerability to future climate change, but this is now a fairly standard approach. It’s\nimportant to critically assess whether such methods are suitable for this species and\nwhether they provide accurate predictions. The genetic offset analysis assumes\nenvironmental adaptation, but in this case, geographic distance is strongly correlated\nwith environmental distance. This makes it difficult to disentangle the effects of\nselection-driven environmental adaptation from neutral genetic drift due to geographic\nseparation. To address this issue, the authors could focus on specific identified adaptive\nvariants and test whether they are associated with environmental variables, geographic\ndifferences, or a combination of both. Differentiating these factors is essential for\nunderstanding the role of environmental adaptation before conducting genetic offset\nanalyses.\nResponse: We thank the reviewer for this thoughtful suggestion. We agree that genetic\nvulnerability analyses are now a widely applied and relatively standard approach. We\nconsidered the approach suitable for B. insignis, a perennial, tree-like fern, and to our\nknowledge this represents the first application of this approach in ferns.\nTo guard against IBD confounding and to ensure that offsets reflect putative\nenvironmental adaptation, we implemented a two-step strategy and have clarified these\ndetails in the revision (Lines 614–617). First, we partitioned variants into neutral and\nadaptive sets. Putatively adaptive SNPs were identified using structure-aware scans\n(LFMM) and redundancy analysis (RDA) with climatic predictors; both methods\naccount for background structure, and we retained only SNPs supported across\napproaches. We then evaluated isolation by environment (IBE) versus isolation by\ndistance (IBD) using Mantel and partial Mantel tests separately for neutral and adaptive\nsets. The adaptive set retained strong environment–genotype associations after\ncontrolling for geographic distance, supporting a role for environment-linked selection\nbeyond spatial drift. Second, we computed genetic offset exclusively on the adaptive\nSNP set, thereby focusing the analysis on loci most plausibly involved in environmental\nresponse. We interpret the resulting offsets as relative genomic vulnerability under the\nclimate scenarios considered and explicitly note caveats relevant to B. insignis (e.g.,\npotential microclimatic buffering, limits to long-distance spore dispersal, and\nunmodeled genotype-by-environment interactions).\nThese clarifications make explicit that (i) the offset analysis is conditioned on evidence\nfor environment-associated adaptation and (ii) the conclusions are framed as\ncomparative risk rather than precise prediction. We thank the reviewer for prompting\nthese improvements in rigor and transparency.\n6. Quantification of Genetic Offset Results: The genetic offset analysis currently\npresents results in the form of a colored plot, which makes it difficult to pinpoint which\npopulations are more vulnerable. I recommend quantifying the offset values (e.g.,\nlocal, forward, and reverse offsets) for each population or genetic cluster and visualizing them\nusing boxplots or similar approaches. This would allow readers to clearly identify\nwhich populations are better adapted and which are more vulnerable to future climate\nchange.\nResponse: We agree that integrating local, forward, and reverse offsets into a single\nmulticolored plot may obscure some of the underlying information. In the revision, we\ntherefore quantified genetic offset at the population level for each metric—local,\nforward, and reverse—and visualized them separately, which are provided in the\nSupplementary Materials (Supplementary Figs. 25-26). The corresponding text has\nbeen revised in the manuscript (Lines 372-390). This revision provides a more complete\npicture and allows readers to more clearly identify which populations are relatively well\nadapted and which may be more vulnerable under future climate scenarios. We\nsincerely appreciate this experienced recommendation, which has improved the clarity\nand interpretability of our results.\nMinor Comments\nOn line 241, the term “core variants” is unclear. How is it defined? Providing a clear\ndefinition would improve clarity.\nResponse: We thank the reviewer for this pointing this out. As shown in Supplementary\nFig. 12, “core variants” refers to the initial set of high-quality SNPs obtained after\nstringent filtering. All downstream datasets (Datasets1–3) were derived from this\nfiltered SNP set. The term was used purely as a technical label without any special\nbiological meaning. To avoid confusion, we have now clarified this definition in the\nMethods section (Lines 535-537).\nIn general, the Results section contains several ambiguous statements, partly because\nthe methods are described later.\nRestructuring or clarifying these sections would make the study easier to follow.\nResponse: We appreciate this helpful comment. We have revised the Results to\neliminate ambiguity arising from deferred methods. Specifically, we (i) define technical\nterms at first mention and provide brief one-sentence method primers, (ii) add explicit\nforward references to the corresponding Methods subsections, and (iii) align the order\nof Results subsections with the analysis workflow. For example, when introducing the\nhigh-confidence SNP set (formerly “core variants”; Lines 240–245), we now define it\nin situ and point to details described in the Methods. Throughout, figure/table callouts\nare synchronized with the text and acronyms are expanded at first use. These changes\nimprove readability without duplicating Methods content.\nResponses to Reviewers’ comments\nReviewer #1 (Remarks to the Author):\nI believe the authors have nicely addressed most of my previous concerns. I found the\nnew synteny comparison with Alsophila very interesting. I do want to point out that\ncomputing “Ks” between LTR repeats is not valid. Ks is synonymous substitution and\nwould only apply to protein-coding genes. How can you calculate Ks for repeats?\nResponse: We appreciate the reviewer’s careful reading. You are absolutely right: Ks\n(synonymous substitutions per synonymous site) applies only to protein-coding genes.\nIn our study, we used LTR_retriever, which applies the Jukes-Cantor (JC69) model to\nestimate the nucleotide substitution rate for noncoding sequences, thereby inferring\nLTR insertion times. The corresponding JC69 formula is , where\n3 4𝐷𝐷\nK is the corrected nucleotide substitution rate rather than Ks𝐾𝐾. W=e− in4aldnv (e1rte−ntl3y) referred\nto this as“Ks”in the text. This change is strictly terminological that does not alter any\nresults or interpretations, and we have revised the wording consistently throughout the\nMethods section (Lines 502-506).\nReviewer #3 (Remarks to the Author):\nI appreciate the efforts of the authors in addressing my previous comments, most of\nwhich have been resolved. Thank you for your efforts. However, I have two minor\ncomments below:\n(1) Although the analysis of genetic load in ferns is limited, it would still be interesting\nto compare the genetic load levels with those of other plant species, such as\nangiosperms. I believe genetic load data are available for more than 20 species, and\nit would be valuable to compare the genetic load of endangered ferns with other\nendangered plant species across phylogenetic positions.\nResponse: We thank the reviewer for this valuable suggestion. Given that genetic\nbackground (e.g., differences in substitution rates), sequencing strategies, data\nprocessing pipelines, and methods for measuring genetic load can vary considerably\namong studies, adopting an appropriate and comparable approach is crucial for\ninterspecific comparisons. To enable fair cross-species comparisons amid variations in\npipelines and ascertainment biases, we employed the widely used proxy π₀/π₄\n(nonsynonymous to synonymous diversity), which reflects the efficacy of purifying\nselection and correlates with genetic load.\nWe compiled π₀/π₄ ratios from >40 plant species spanning ferns and seed plants, and\nadded a new Supplementary Figure summarizing these values (Supplementary Fig. 21).\nThe results show that the endangered ferns (Brainea insignis and Alsophila spp.) exhibit\nsignificantly elevated π₀/π₄ relative to most seed plants, with the exception of Acer\nyangbiense, a critically endangered species with extremely small effective population\nsize (Nₑ). These patterns are consistent with the expectation that small effective\npopulation sizes weaken purifying selection, thereby increasing the burden of\ndeleterious polymorphisms.\nWe acknowledge that future studies incorporating more genomic data may allow for\nfurther updates or refinements of these results. Nevertheless, we believe that presenting\nthis comparative analysis in the supplementary materials provides valuable contextual\ninformation and serves as a useful reference for understanding the genetic load\ncharacteristics of endangered ferns relative to other plant groups.\nSupplementary Fig. 21 | Summary of π₀/π₄ ratios across fern and seed plant species.\n(2) The forward and reverse differences across different populations shown in\nSupplementary Figure 25 are difficult to distinguish, likely because the local offset\nvalues are too high. I suggest that instead of using the same value range for all the\nplots, each figure should use its own value range for plotting the offset. This\nadjustment would help to clearly visualize the landscape distribution of offset\nvalues and make it easier to differentiate between forward and reverse offsets.\nResponse: Thank you for the suggestion. We have adjusted the plotting ranges so that\nthe local, forward, and reverse offset maps now each use their own respective value\nranges, instead of the previous approach that applied a range based solely on the local\noffset values (Supplementary Figure 26). This adjustment has indeed significantly\nenhanced the visual distinguishability between forward and reverse offsets, as\nillustrated below. Importantly, these changes do not affect the original results; they only\nimprove the presentation of the offset landscapes.\nOriginal Supplementary Fig. 25\nRevised Supplementary Fig. 26\nRevision #4\nResponses to Reviewers’ comments\nReviewer #3 (Remarks to the Author):\nI appreciate the authors for addressing my concerns. I have no further comments and\nlook forward to seeing this work published soon.\nResponse: We sincerely thank Reviewer #3 for the constructive feedback in the review\nprocess. Your comments have greatly improved our manuscript and aligned it better\nwith the journal’s standards. We appreciate your support."}, {"round": 3, "reviewer_comments": {"reviewer_3": "I appreciate the authors for addressing my concerns. I have no further comments and look forward to seeing this work published soon."}, "author_response": "Response: We sincerely thank Reviewer #3 for the constructive feedback in the review process. Your comments have greatly improved our manuscript and aligned it better with the journal’s standards. We appreciate your support."}]

Decision letter