Abstract

Background

   Leafdevelopment represents a crucial stage in the plant life cycle,
   involving complex morphogenetic and physiological processes governed by
   evolving molecular mechanisms and metabolite profiles. The growth and
   maturation of Angiopteris fokiensis Hieron, a species used in
   traditional Chinese medicine, are characterized by fluctuating
   metabolite accumulation patterns regulated by largely unknown molecular
   pathways.

Results

   Touncover these pathways, we employed next-generation sequencing to
   construct the A. fokiensis leaf transcriptome at two distinct
   developmental stages, allowing for a comprehensive analysis of gene
   expression dynamics while emphasizing the identification of genes that
   regulate leaf development and metabolite synthesis. The de novo
   assembly of high-quality sequencing reads generated 117,627 unigenes
   averaging 1,308 base pairs in length. FPKM analysis uncovered
   significant transcriptomic alterations during leaf development.
   Additionally, non-targeted metabolomics identified 1,494 distinct
   analytes, with lipids representing the most abundant metabolite class
   in both A. fokiensis samples. In the 'phenylalanine, tyrosine and
   tryptophan biosynthesis' pathway, two downregulated arogenate
   dehydrogenase (NADP+) genes (Unigene23378-S4 and Unigene47537-S2) in
   Stage1 correlated with reduced L-tyrosine levels. In the 'galactose
   metabolism' pathway, the upregulation of three beta-galactosidase genes
   (Unigene43641-S6, Unigene43648-S6, Unigene47074-S1) and the
   downregulation of one (Unigene28294-S2) corresponded to decreased
   alpha-lactose levels.

Conclusions

   This study provides an in-depth examination of the dynamic
   transcriptomic and metabolomic changes occurring during A. fokiensis
   leaf development, revealing key regulatory networks and enhancing the
   annotation of theA. fokiensis genome. These findings lay a crucial
   groundwork for future research on this medicinal plant.

Supplementary Information

   The online version contains supplementary material available at
   10.1186/s12863-025-01366-7.

   Keywords: Angiopteris fokiensis, Leaf, Development, Transcriptome,
   Metabolite

Background

   Leaf development is a vital stage in a plant’s life cycle, encompassing
   a series of complex and dynamic processes governed by numerous genes
   that affect leaf morphology, photosynthetic capacity, and metabolite
   accumulation. Next-generation sequencing (NGS) technologies, such as
   Illumina High-Seq, have transformed genomic research by providing
   high-throughput sequencing at relatively lower costs than traditional
   methods. These platforms are particularly useful for applications like
   de novo sequencing, genome resequencing, and transcriptome analysis
   [[32]1], especially in non-model organisms lacking a reference genome,
   such as Angiopteris fokiensis Hieron.

   NGS data offer valuable insights into the molecular mechanisms
   regulating gene expression [[33]2], providing rich de novo genomic or
   transcriptomic resources for gene discovery, molecular marker
   development, gene localization, gene expression studies, and
   comparative genomics. For example, Illumina paired-end sequencing has
   been employed to investigate dynamic transcriptome changes during leaf
   growth and development in several plant species, including Apium
   graveolens L [[34]3]., Eucommia ulmoides [[35]4], Bergera koenigii L
   [[36]5]., Panicum virgatum L [[37]6]., Osmanthus fragrans [[38]7],
   Brassica campestris L [[39]8]., and Epimedium pubescens [[40]9].

   Metabolite profiling offers an additional approach for investigating
   plant developmental processes by revealing how fluctuations in
   metabolite levels shape phenotypes and capturing changes occurring
   during leaf development [[41]4, [42]9, [43]10]. This approach enables
   the identification of metabolites by linking them to metabolic
   networks, functions, and pathways [[44]11]. Integrating metabolomics
   with transcriptomics has further advanced our understanding of
   biosynthetic pathways involved in producing essential metabolic
   compounds [[45]12, [46]13], particularly in plants that generate
   bioactive compounds with medicinal properties [[47]14, [48]15].

   Similarly, transcriptome and metabolome analyses in ferns have revealed
   distinctive metabolic profiles. Studies have included species such as
   Hymenoglossum cruentum, Hymenophyllum dentatum [[49]16], Cyclosorus
   parasiticus [[50]17], Cyathea delgadii [[51]18], Drynaria roosii
   [[52]19], cloud forest ferns [[53]20], and tree ferns Angioptera
   spinulosa and Angioptera metteniana [[54]21]. Metabolomics offers
   crucial insights into metabolite concentrations and changes in
   metabolic pathways during development. Angiopteris fokiensis, as a
   prized tropical tree fern within the Angiopteridaceae family, is also
   economically important. However, despite advancements in plant
   metabolomics and transcriptomics, the specific biochemical and
   molecular mechanisms regulating A. fokiensis leaf growth and
   development remain poorly understood due to limited data on
   transcriptional changes across its developmental stages.

   This study aims to fill this knowledge gap by utilizing Illumina
   paired-end sequencing to develop a comprehensive leaf transcriptome for
   A. fokiensis at two distinct developmental stages. In parallel,
   metabolome profiling was conducted to categorize and quantify leaf
   metabolites and monitor changes in their accumulation during leaf
   growth and development. Our results provide new insights into the
   molecular processes governing metabolite biosynthesis and regulation in
   A. fokiensis leaves, highlighting the power of an integrated approach
   in advancing our understanding of plant development.

Methods

Plant materials

   A. fokiensis plants were obtained from the Medicine medicinal plant
   resource nursery (Guizhou University of Traditional Chinese) in March
   2024. Transcriptome analysis was performed on representative tissues
   collected at two distinct growth stages, designated as Stage1 and
   Stage2. Both immature and mature leaf tissue (differentiated by size,
   with immature leaves ranging from 4 cm, and mature leaflets from 12 cm)
   were collected from the A. fokiensis. Leaf samples were immediately
   stored in liquid nitrogen at − 190 °C. Figure [55]1 shows both immature
   and mature leaves from the A. fokiensis used in this study.
   Transcriptome analysis for each stage involved three biological
   replicates. Metabolome analysis for each stage involved six biological
   replicates.

Fig. 1.

   [56]Fig. 1
   [57]Open in a new tab

   Morphological characterization of A. fokiensis leaves at different
   growth stages

RNA extraction and illumina sequencing

   Total RNA was isolated from leaves using the plant total RNA extraction
   kit (TIANGEN). Quality control (QC) analysis was next performed to
   evaluate RNA preparation quality based on RNA concentration, RNA
   integrity number (RIN), 28 S/18S ratio, and fragment length
   distribution using an Agilent 2100 Bioanalyzer (Agilent RNA 6000 Nano
   Kit). RNA purity was assessed using a NanoDrop™ spectrophotometer. The
   RNA samples were sent to the Beijing Genomics Institute (BGI) for
   further purification, library construction, and RNA-Seq analysis.

   Raw sequencing data were processed with SOAPnuke (v1.6.5) to remove
   low-quality reads, and the remaining clean reads were saved in FASTQ
   format [[58]22]. These reads were assembled into contigs using Trinity
   (v2.13.2), and the assembly quality was assessed using the Benchmarking
   Universal Single-Copy Orthologs (BUSCO) tool [[59]23]. Gene expression
   levels were quantified using RNA-Seq by Expectation-Maximization
   software (RSEM, v1.3.1) [[60]24]. Protein-coding sequences within
   unigenes were identified using TransDecoder (v5.5.0) [[61]25] and
   aligned to the SwissProt database using the Basic Local Alignment
   Search Tool (BLAST). Hmmscan was employed to predict Pfam protein
   domains and coding regions.

   Functional annotation of unigenes, including transcription factor (TF)
   functions, was performed using the Gene Ontology (GO) and Kyoto
   Encyclopedia of Genes and Genomes (KEGG) databases. Differentially
   expressed genes (DEGs) were analyzed with DESeq2 (v1.4.5) [[62]26],
   using a significance cutoff of Q ≤ 0.05 (or FDR ≤ 0.001). Visualization
   of DEG patterns was performed by creating a heatmap of DEG clusters
   using the pheatmap function. Post-analysis, DEGs were functionally
   annotated using the hypergeometric distribution method (phyper function
   in R) for KEGG pathway enrichment
   ([63]https://en.wikipedia.org/wiki/Hypergeometric_distribution) and the
   TermFinder package for GO enrichment
   ([64]https://metacpan.org/pod/GO::TermFinder). DEGs with significant
   enrichment were identified based on a Q value threshold of ≤ 0.05.

Metabolic profiling

   Samples (six biological replicates per stage) were subjected to
   non-targeted metabolome analysis. Equipment for sample preparation
   included a vortex mixer (QL-901, Qilinbeier Instrument Manufacturing
   Co., Ltd.), low-temperature high-speed centrifuge (Eppendorf 5430),
   refrigerated vacuum concentrator (Maxi Vac Beta, GENE Co.), Milli-Q
   water purification system (Integral Millipore Corp., USA), and a
   weaving tissue grinder (JXFSTPRP). Reagents included LCMS grade
   methanol (A454-4) and acetonitrile (A998-4) (Thermo Fisher Scientific,
   USA), ammonium formate (17843-250G, Honeywell Fluka, USA), and formic
   acid (50144-50 ml, DIMKA, USA). Internal standards included d3-Leucine,
   13C9-Phenylalanine, 13C3-Progesterone, and d5-Tryptophan.

   Metabolite extraction was initiated by placing each 50-µg sample in a
   separate 1.5-mL Eppendorf tube, followed by the addition of 800 µL of
   precooled extraction solution (methanol: H[2]O = 7:3, v/v) and 20 µL of
   internal standard 1 (IS1). Samples were homogenized with a weaving
   tissue grinder for 10 min at 50 Hz, then sonicated for 30-min in a 4 °C
   water bath. Tubes were incubated upright without agitation at −20 °C
   for 1 h then extracts were centrifuged at 4 °C for 15 min at 14,000 rpm
   and supernatants (600 µL) were filtered through 0.22-µm filters. To
   ensure LC-MS repeatability and stability, 20-µL volumes of filtered
   supernatants were pooled to create quality control (QC) samples. Both
   individual and pooled QC samples were then transferred to 1.5-mL vials
   for LC-MS analysis.

   LC-MS analysis was conducted using an ultra-high-performance liquid
   chromatography (UPLC) system (I-Class Plus, Waters, USA) equipped with
   a Hypersil GOLD aQ column (1.9 μm, 100 mm × 2.1 mm, Thermo Fisher
   Scientific) and a Q Exactive high-resolution mass spectrometer (Thermo
   Fisher Scientific). Each sample underwent chromatographic separation
   with an injection volume of 5 µL, a column temperature of 40 °C, and
   mobile phases consisting of 0.1% formic acid in water (Phase A) and
   0.1% formic acid in acetonitrile (Phase B). The mobile phases were
   delivered at a flow rate of 0.3 mL/min under the following gradient
   conditions: 5% B over 0.0–2.0 min, 5–95% B over 2.0–22.0 min, hold at
   95% B over 22.0–27.0 min, and washing with 95% B over 27.1–30 min.

   Mass spectrometry (MS) data were acquired from primary and secondary
   scans ranging from 125 to 1500 m/z for positive ions and 100–1500 m/z
   for negative ions at a resolution of 70,000. The automatic gain control
   (AGC) target was set to 1e6 with a maximum ion injection time of 100
   ms. For MS/MS fragmentation, the top three precursors were selected at
   a resolution of 30,000. The AGC target was set to 2e5 with a maximum
   ion injection time of 50 ms using stepped normalized collision energy
   settings of 20, 40, and 60 eV. The electrospray ionization (ESI)
   parameter settings included a sheath gas flow rate of 40, auxiliary gas
   flow rate of 10, positive-ion mode spray voltage (|KV|) of 3.80,
   negative-ion mode spray voltage (|KV|) of 3.20, capillary temperature
   of 320 °C, and auxiliary gas heater temperature of 350 °C.

   MS data were processed using Compound Discoverer 3.3 (Thermo Fisher
   Scientific) and analyzed using the BGI metabolome database (BMDB),
   ChemSpider online database, and the mzCloud database, resulting in the
   generation of a data matrix encompassing metabolite peak areas and
   identification results. Further data processing was performed using
   Compound Discoverer (v.3.3) with the following settings: parent ion
   mass deviation < 5 ppm, fragment ion mass deviation < 10 ppm, and
   retention time deviation < 0.2 min
   ([65]https://mycompounddiscoverer.com/).

   Data normalization was performed using probabilistic quotient
   normalization (PQN) [[66]27], with QC-based robust locally estimated
   scatterplot smoothing (LOESS) signal correction (QC-RLSC) [[67]28]
   employed to adjust for batch effects. Metabolites with a coefficient of
   variation exceeding 30% in the QC samples were excluded. After log2
   transformation, a partial least squares-discriminant analysis (PLS-DA)
   model was constructed and employed for comparative analysis of the two
   sample groups performed using Pareto scaling and 7-fold
   cross-validation. Orthogonal projections to latent structures
   discriminant analysis (OPLS-DA) was also applied to decompose the data
   matrix into components related and unrelated to the response variable,
   enhancing model simplicity, predictive accuracy, and explanatory power.

Real-time quantitative PCR analysis

   Real-time quantitative PCR was conducted to validate transcript levels
   of genes identified by RNA-Seq. The analysis was performed using AceQ^®
   qPCR SYBR Green Master Mix (Vazyme) on a RealTime PCR System (Roche
   LightCycler 480 II, 384). RNA samples from A. fokiensis, previously
   used for RNA-Seq, were reverse transcribed into cDNA, and each qPCR
   reaction was performed in three technical replicates. The specificity
   of the amplification was confirmed through melting curve analysis and
   gel electrophoresis of the final products. Cycle threshold (CT) values
   for each gene were normalized against a reference gene, and relative
   gene expression levels were calculated using the ΔΔC[T] method as
   outlined in previous studies [[68]29].

Results

DEG analysis in A. fokiensis leaves during two developmental stages

   To investigate the molecular mechanisms underlying A. fokiensis leaf
   growth and development, a comprehensive transcriptome analysis was
   performed using an RNA-Seq Analyzer II system. This analysis aimed to
   identify DEGs between two developmental stages of A. fokiensis leaves
   (Fig. [69]1), utilizing six cDNA Libraries generated from total RNA
   samples. Each sample produced over 45.12 million raw reads, resulting
   in approximately 43.87 million high-quality (clean) reads per sample
   (Table [70]1). Data quality, confirmed based on a Q30 score of ≥ 92.8%
   for each sample, demonstrated data reliability for further analysis.

Table 1.

   Summary of RNA-Seq datasets for the six libraries
   Samples Raw Reads (M) Clean Reads (M) Clean Bases (Gb) Clean Reads Q20
   (%) Clean Reads Q30 (%)
   Stage1-1 46.08 44.26 6.64 97.67 92.8
   Stage1-2 45.76 44.13 6.62 97.72 93.02
   Stage1-3 45.12 43.87 6.58 97.72 92.96
   Stage2-1 47.68 44.23 6.63 97.88 93.59
   Stage2-2 49.04 44.09 6.61 98.02 94.03
   Stage2-3 49.6 44.22 6.63 97.93 93.78
   [71]Open in a new tab

   High-quality reads from two developmental stages were combined to
   generate a comprehensive leaf transcriptome dataset. Short reads were
   de novo assembled using Trinity software, producing 117,627 unigenes
   (Fig. S1) averaging 1,308 bp in length. Open reading frames (ORFs) were
   identified using TransDecoder, yielding a collection of nucleotide and
   protein sequences suitable for gene cloning, phylogenetic analysis, and
   functional verification. Among these unigenes, 24,705 contained
   high-integrity ORFs with complete coding sequences (CDSs) from
   initiation to stop codons (Fig. S2). The transcriptomic profiles of the
   different developmental stages showed distinct patterns, with
   consistent results obtained across biological replicates.

   Gene expression changes based during the two developmental stages were
   quantified as fragments per kilobase of transcript per million mapped
   reads (FPKM). To validate the RNA-Seq findings, a subset of genes was
   selected for quantitative real-time (qRT)-PCR analysis, which was
   conducted using seven primer pairs, including one for the reference
   gene encoding A. fokiensis actin (Table S1). The correlation between
   RNA-Seq and RT-PCR data showed strong agreement in gene expression
   trends (Supplemental Fig. S3), with a correlation heat map further
   confirming high consistency among the three biological replicates (Fig.
   S4).

   DEG analysis revealed significant transcriptomic changes during leaf
   development in A. fokiensis. A total of 992 genes were up-regulated and
   935 were down-regulated in Stage2 compared to Stage1, with FPKM values
   provided in Table S2. This analysis highlighted distinct gene
   expression patterns between the two stages, suggesting that a complex
   regulatory mechanism modulates gene expression to regulate signal
   transduction during leaf development. Additionally, GO (Fig. S5) and
   KEGG (Fig. [72]2 and Table S3) enrichment analyses of DEGs identified
   key biological processes and pathways associated with the different
   developmental stages of A. fokiensis leaves, including “metabolic
   pathways,” “biosynthesis of secondary metabolites,” and “carbon
   metabolism.”

Fig. 2.

   [73]Fig. 2
   [74]Open in a new tab

   KEGG pathway enrichment analysis of all genes expressed in two
   different stages of A. fokiensis leaf development

Metabolic differences observed between two A. fokiensis leaf developmental
stages

   Principal component analysis (PCA) of the metabolomic data showed a
   clear separation between Stage1 and Stage2 leaves along PC1 (53.41%)
   and PC2 (14.49%) (Fig. [75]3A). Ultimately, a total of 1,494
   metabolites were identified and grouped into thirty distinct categories
   (Fig. S6). The most abundant compounds included “lipids” (104),
   followed by “others” (36), “terpenoids” (35), “benzene and derivatives”
   (27), “carbohydrates” (22), “alkaloids” (21), and “organic acids” (16).
   Additionally, the metabolites were classified into 18 categories
   according to KEGG analysis findings (Fig. S7).

Fig. 3.

   [76]Fig. 3
   [77]Open in a new tab

   Characterization of metabolite profiles in A. fokiensis leaves at two
   developmental stages. A Principal Component Analysis (PCA) score plot
   illustrating the separation of metabolite profiles between Stage 1 and
   Stage 2 leaves. B Number of significantly changed metabolites (SCMs)
   between Stage 1 and Stage 2 leaves. C Heatmap and cluster analysis of
   metabolite profiles at the metabolome level, comparing different stages
   of leaf development. D KEGG pathway annotation of SCMs identified
   between the two stages of A. fokiensis leaves

   Significantly changed metabolites (SCMs) between the two stages were
   identified by analyzing six biological replicates of leaf samples
   (Fig. [78]3B, C and Table S4). A total of 578 metabolites showed
   significant differences in accumulation between the two stages, with
   316 increased and 262 decreased levels in Stage2 compared to Stage1.
   KEGG pathway enrichment analysis of SCMs revealed pathways related to
   “biosynthesis of cofactors”, “carbon metabolism”, and “glyoxylate and
   dicarboxylate metabolism”, among others (Fig. [79]3D).

Integrative transcriptome and metabolome analysis

   To integrate metabolomics with transcriptomics data, we identified DEGs
   and metabolites involved in key biological pathways. This approach
   aimed to identify functional genes or pathways jointly enriched in
   metabolites, thereby defining pathways linked to crucial biological
   phenotypes. The association analysis used quantitative values of
   differential metrics for a combined assessment.

   Integrative analysis revealed DEGs and metabolites associated with key
   biological pathways. By mapping DEGs and SCMs to KEGG pathways, a
   comprehensive network was constructed, uncovering functional
   relationships between gene expression and metabolite accumulation. This
   analysis identified pathways Linked to crucial biological phenotypes,
   incorporating quantitative metrics for a combined assessment. DEGs were
   enriched in 277 pathways, and SCMs in 66 pathways, with several
   intersecting in critical pathways.

   Pathway bubble plots generated from enrichment analysis identified
   significantly enriched pathways (Fig. [80]4), including ‘ascorbate and
   aldarate metabolism’ (map00053) (Fig. S8), ‘arginine biosynthesis’
   (map00220) (Fig. S9), ‘alanine, aspartate and glutamate metabolism’
   (map00250) (Fig. S10), ‘arginine and proline metabolism’ (map00330)
   (Fig. S11), and ‘alpha-linolenic acid metabolism’ (map00592) (Fig.
   S12). For example, in the ascorbate and aldarate metabolism pathway,
   133 genes encoding relevant enzymes were identified, with 12 DEGs and
   one SCM showing strong correlations.

Fig. 4.

   [81]Fig. 4
   [82]Open in a new tab

   Pathway enrichment analysis bubble plot. X-axis enrichment factor
   (richfactor), the larger the value illustrates the greater the
   proportion of differential metabolites annotated to that pathway and
   the differential metrics being associated with. The circles represent
   the indicator pathway of the pathway being associated with omics, the
   triangles represent metabolic pathways, the graph size represents the
   number of differential metabolites annotated to that pathway versus the
   indicator of the difference in the pathway being associated with omics,
   and the graph color indicates pathway significance

   Predictions of gene-regulated metabolite changes were also made. For
   example, in the ‘phenylalanine, tyrosine and tryptophan biosynthesis’
   pathway (map00400), two downregulated arogenate dehydrogenase (NADP+)
   genes (Unigene23378-S4 and Unigene47537-S2) in Stage1 correlated with
   reduced L-tyrosine levels, implicating their role in ascorbate
   biosynthesis (Fig. [83]5). Similarly, in the ‘galactose metabolism’
   pathway (map00052), the upregulation of three beta-galactosidase genes
   (Unigene43641-S6, Unigene43648-S6, Unigene47074-S1) and the
   downregulation of one (Unigene28294-S2) corresponded to decreased
   alpha-lactose levels.

Fig. 5.

   [84]Fig. 5
   [85]Open in a new tab

   Regulation of metabolites by genes. Red and green boxes represent
   upregulated and downregulated genes, respectively. Red and green dots
   indicate metabolites with increased and decreased accumulation,
   respectively

Discussion

   A. fokiensis is widely distributed in southern China, where it holds
   significant medicinal, ecological, and economic value. Beyond its
   traditional medicinal applications, it plays a role in landscaping and
   soil and water conservation. However, over-harvesting and habitat
   destruction have led to a dramatic decline in wild populations. Despite
   its medicinal potential and market demand, the molecular mechanisms
   regulating metabolite biosynthesis during its leaf growth and
   development remain largely unexplored. This study addresses this gap by
   integrating transcriptomic and metabolomic analyses to elucidate the
   regulatory networks and metabolite dynamics of A. fokiensis leaves at
   two developmental stages. In the absence of a published genome for A.
   fokiensis, the annotation of nearly 32.02% of unigenes against known
   genes or proteins across multiple databases represents a notable
   achievement. High correlation coefficients among biological replicates,
   coupled with confirmatory qPCR results, underscore the reliability of
   the transcriptomic data. Differential expression analysis revealed
   substantial transcriptomic changes, with DEGs significantly enriched in
   KEGG pathways such as ‘biosynthesis of secondary metabolites’,
   ‘metabolic pathways’, and ‘carbon metabolism’. Similar enrichment
   patterns have been observed in other plant species, such as B.
   campestris L. (Shi et al., 2023), where DEGs were linked to processes
   like ‘plant hormone signal transduction’, ‘starch and sucrose
   metabolism’, and ‘circadian rhythm’, and in E. pubescens (Xu et al.,
   2023), which showed significant enrichment in pathways such as ‘plant
   hormone signal transduction’, ‘MAPK signaling pathway – plant’, and
   ‘RNA transport’. In this study, transcription factors (TFs) emerged as
   critical regulators of these pathways (Table [86]S2) (Fig. S13). For
   instance, genes such as Unigene46787-S1 (Phoenix dactylifera
   vacuolar-sorting receptor 3, WRKY) and Unigene58608-S4 (Camptotheca
   acuminata isolate Cac002 AP2/ERF TF, AP2) exhibited higher expression
   in Stage2 compared to Stage1, while Unigene20299-S6 (Theobroma cacao
   telomere repeat-binding protein 5, MYB) and Unigene38653-S4 (Gossypium
   raimondii expansin-A15-like, bHLH) showed lower expression in Stage2.
   These results highlight TFs with distinct expression profiles during
   different developmental stages of A. fokiensis leaves, underscoring
   their stage-specific roles. However, further studies are required to
   confirm the functions of these TF genes.

   The widespread differential regulation of genes and metabolites across
   various metabolic pathways underscores the complexity of A. fokiensis
   leaf development, characterised by multiple transitions governed by a
   network of interacting genes and signaling pathways. Metabolic
   profiling identified 1,494 metabolites across two developmental stages,
   classified into 30 categories. Significant metabolite classes included
   ‘amines’ (e.g., tyramine, oleoyl ethanolamide), ‘cofactors’ (e.g.,
   5,6,7,8-tetrahydromethanopterin), ‘fatty acyls’ (e.g., anandamide,
   palmitoyl ethanolamide), ‘indazoles’ (e.g., granisetron), ‘lignans’
   (e.g., podorhizol beta-D-glucoside), and ‘organic acids’ (e.g.,
   dicyclomine, 2-butoxyacetic acid, citric acid, 2-oxoglutaric acid). The
   significant accumulation of these metabolites in growing leaves
   highlights their role in supporting active growth and development
   (Table S4). In comparison, a study on E. ulmoides, performed via
   UPLC-MS, identified only 515 metabolites, including 127 flavonoids, 46
   organic acids, 44 amino acid derivatives, 8 isoflavones, 9
   phenolamides, and 16 vitamins (Li et al., 2019). Similarly, a
   metabolomic analysis in rice detected only 510 and 512 metabolites in
   two rice varieties [[87]24]. Among the 578 significantly changed
   metabolites (SCMs) identified between developmental stages, 316 showed
   increased accumulation in Stage2 compared to Stage1, whereas 262 showed
   decreased accumulation. These changes Likely influence leaf
   morphogenesis, photosynthetic capacity, and metabolic activity by
   modulating biosynthetic and degradative pathways. Notable SCMs included
   chlortetracycline, kaempferol
   3-neohesperidoside-7-(2’’-p-coumaryllaminaribioside), and dermorphin,
   which exhibited significant variations. Pathway enrichment analysis
   linked these SCMs to key pathways such as ‘biosynthesis of cofactors’,
   ‘carbon metabolism’, and ‘glyoxylate and dicarboxylate metabolism’.
   Integrative transcriptomic and metabolomic analysis provided deeper
   insights into the regulatory networks governing A. fokiensis leaf
   development. For example, in the ‘ascorbate and aldarate metabolism’
   pathway, 12 DEGs and one SCM were identified, highlighting coordinated
   regulation during development. Similarly, in the ‘phenylalanine,
   tyrosine, and tryptophan biosynthesis’ pathway, the downregulation of
   two arogenate dehydrogenase genes corresponded with decreased
   L-tyrosine levels in Stage1, suggesting developmental regulation of
   ascorbate biosynthesis. This study represents a significant step
   towards understanding the molecular basis of A. fokiensis leaf
   development and metabolite biosynthesis. The findings lay the
   groundwork for functional genomics research and may inform conservation
   and sustainable utilisation strategies for this valuable medicinal
   plant. However, limitations such as the small sample size and lack of a
   reference genome highlight the need for further investigations to
   validate the identified regulatory networks and pathways. Despite these
   limitations, the results demonstrate the potential of genome-wide
   approaches to unravel the molecular mechanisms underlying A. fokiensis
   development and evolution.

   Leaf development is a vital stage in a fern’s life cycle, encompassing
   a series of complex and dynamic processes governed by numerous genes
   and metabolite. The physical protection of young leaves is adapted to
   their vulnerable growth stage, whereas the chemical defense and
   reproductive functions of mature leaves ensure population persistence.
   Investigating these differences may provide insights into fern
   evolution and medicinal development.

Conclusions

   This study offers dynamic models of transcriptomic and metabolite
   changes, uncovering key regulatory networks driving A. fokiensis leaf
   growth and development. By enhancing gene annotation for A. fokiensis,
   it establishes a robust foundation for future research into the
   molecular mechanisms underlying its development. These insights hold
   significant promise for advancing its applications in traditional
   medicine and promoting sustainable utilisation.

Supplementary Information

   [88]Supplementary Material 1.^ (1.8MB, xls)
   [89]Supplementary Material 2.^ (1.6MB, pdf)

Acknowledgements