The genus Phoenix (Arecaceae) comprises economically, eco logically, and evolutionarily significant palm species that exhibit remarkable diversity in growth form, ecological tolerance, and geographic distribution. While Phoenix dactylifera (date palm) has been extensively studied at the genomic [2] and transcrip tomic levels [42], comparatively little is known about the mo lecular basis underlying functional divergence among wild and semi-wild Phoenix species. Understanding transcriptional varia t ion among closely related species provides valuable insights into the genetic mechanisms associated with ecological adapta t ion, physiological specialization, and trait evolution [40].

Comparative transcriptome analysis using RNA sequencing (RNA-seq) has emerged as a powerful approach to investigate gene expression divergence across species, particularly in non model plants where complete genome sequences may be un available [8, 30]. Reference-guided alignment strategies using high-quality genomes from related species have been success fully employed to study transcriptional variation in wild rela t ives of crop plants, enabling robust differential expression and functional annotation analyses [4,16]. Such approaches have been widely applied to study leaf transcriptomes, given the central role of leaves in photosynthesis, stress perception, and metabolic regulation [11,39].

Leaves are particularly informative tissues for comparative studies because they integrate developmental, physiological, and environmental signals. Comparative leaf transcriptomics has been used to uncover adaptive gene expression patterns related to photosynthesis efficiency, secondary metabolism, water-use efficiency, and abiotic stress tolerance across diverse plant lineages [9,26]. Differences in leaf gene expression be tween species often reflect ecological specialization and evo lutionary divergence, even when overall genomic similarity re mains high [34].

P. theophrasti and P. pusilla represent two distinct evolution ary lineages within the genus, differing in geographic distribu t ion, ecological preferences, and morphological traits. While P. theophrasti is restricted to the eastern Mediterranean region, P. pusilla is adapted to semi-arid and coastal environments of South Asia, where it persists under contrasting climatic and edaphic conditions. Despite their close phylogenetic relation ship, the molecular mechanisms underlying their ecological dif ferentiation remain largely unexplored [6].

In this study, we performed a comparative transcriptome analysis of mature leaf tissues from P. theophrasti and P. pu silla using a reference-guided RNA-seq framework based on the P. dactylifera genome. Differential gene expression analysis, Gene Ontology (GO) enrichment, and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were conducted to identify species-specific transcriptional signatures. Further more, we curated differentially expressed genes associated with key functional traits to highlight biologically meaningful divergence between the two species. Our results provide novel insights into transcriptional differentiation within wild Phoenix species and contribute to a broader understanding of functional evolution in palms.

Plant material and RNA isolation

Mature leaf tissues were collected from plants of P. theophrasti and P. pusilla growing under field conditions at the USDA-Subtropical Horticulture Research Station (SHRS), Miami, Florida (latitude N25°38’33.76” and longitude 80°17’37.86). For each species, leaves from three individuals were harvested to capture representative transcriptional profiles and pooled be fore RNA extraction. Total RNA was extracted from plant tissue using the RNeasy Plant Mini Kit (Qiagen, Hilden, Germany) in accordance with the manufacturer’s recommended protocol to ensure high-quality RNA recovery. The concentration and purity of the isolated RNA were initially evaluated using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, MA, USA), which measures absorbance ratios to assess contamination. Further assessment of RNA integrity was performed using an Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA), pro viding an RNA Integrity Number (RIN) as a standardized mea sure of RNA quality. Only samples exhibiting a RIN value greater than 7.0, indicative of intact and high-quality RNA, were select ed for subsequent library preparation and sequencing steps.

RNA sequencing and quality control

RNA sequencing libraries were constructed with the Kapa Hyper Stranded mRNA library kit (Roche) following the manu facturer’s guidelines to preserve strand specificity. The process involved enrichment of poly(A)-tailed mRNA, fragmentation into smaller pieces, and reverse transcription into complemen tary DNA (cDNA). After cDNA synthesis, sequencing adapters were ligated to the fragments, and the libraries were amplified to generate sufficient material for sequencing. High-quality li braries were then subjected to next-generation sequencing (NGS) at the University of Illinois Roy J. Carver Biotechnology Center. Sequencing was performed on the Illumina NovaSeq 6000 platform (Illumina, CA, USA) using the NovaSeq SP re agent kit, producing paired-end reads of 150 base pairs (bp), which provide comprehensive coverage and improved accuracy for downstream transcriptomic analyses. Raw sequencing reads were subjected to quality assessment using FastQC [3]. Adapter sequences and low-quality bases were removed using Trimmo matic, retaining reads with high base quality for downstream analyses [5]. All sequence reads were deposited in the public NCBI Sequence Read Archive database under BioSample acces sions - PRJNA1400048.

Reference-guided read alignment

Cleaned reads from both species were aligned to the P. dac tylifera reference genome using HISAT2, a splice-aware aligner optimized for RNA-seq data [23]. The reference-guided align ment strategy was employed due to the close phylogenetic re lationship among Phoenix species and the availability of a high quality P. dactylifera genome [2]. Default alignment parameters were used unless otherwise specified. Alignment statistics were examined to ensure mapping quality and consistency across samples.

Transcript quantification and gene expression estimation

Aligned reads were assigned to annotated gene features us ing featureCounts from the Subread package [27]. Raw gene level read counts were generated for each species and used as input for differential expression analysis. Genes with extremely low read counts were filtered before downstream statistical analyses to reduce noise.

Differential gene expression analysis

Differential gene expression analysis between P. theophrasti and P. pusilla leaf transcriptomes was performed using DESeq2 [28]. Read counts were normalized using DESeq2’s median of-ratios method. Differentially expressed genes (DEGs) were identified based on log2 fold change (Log2 FC) absolute value of ≥1 and adjusted p-values (0.05) calculated using the Ben jamini–Hochberg False Discovery Rate (FDR) correction. Given the pooled nature of the samples, DEGs were interpreted as species-level transcriptional differences rather than population level variation.

Functional annotation and gene ontology enrichment analysis

Functional annotation of expressed genes was derived from the P. dactylifera reference genome annotations. Gene Ontol ogy (GO) enrichment analysis of DEGs was performed using the clusterProfiler R package [45]. Enrichment analyses were con ducted separately for Biological Process, Molecular Function, and Cellular Component categories. GO terms with adjusted p-values below the selected significance threshold were consid ered significantly enriched.

KEGG pathway enrichment analysis

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was conducted using clusterProfiler [45] in conjunction with KEGG annotations [21]. Differentially ex pressed genes were mapped to KEGG pathways to identify bio logical pathways showing significant enrichment between the two Phoenix species. Enriched pathways were interpreted in the context of leaf physiology, metabolism, and species-specific functional traits.

Trait-associated gene identification

To identify biologically relevant transcriptional differences, DEGs were manually curated into functional trait categories based on annotation descriptions, GO terms, and KEGG path way associations. Trait-associated gene sets were selected to highlight processes relevant to leaf function, metabolism, and stress-related pathways. These curated gene lists were used for downstream visualization and interpretation.

Data visualization and statistical analysis

All statistical analyses and visualizations were performed in R. MA plots, volcano plots, heatmaps, and enrichment plots were generated using a combination of DESeq2, clusterProfiler, and ggplot2 [41]. Heatmaps were produced using Pheatmap v1.0.12 [24] with the scaled normalized expression values to visualize relative gene expression patterns between P. theo phrasti and P. pusilla.

Sequencing output and reference-guided alignment

RNA-seq libraries generated from mature leaf tissues of P. theophrasti and P. pusilla produced high-quality paired-end se quencing data with consistent depth across samples. Following quality filtering, reads were aligned to the P. dactylifera refer ence genome using HISAT2. Overall alignment rates were high and comparable between species, ranging from approximately 84% to 85% across all libraries (Table 1). The number of ex pressed genes detected was also similar between species, with approximately 26,000–27,000 genes expressed per library.

The consistently high alignment efficiency across both spe cies indicates that the P. dactylifera reference genome provides adequate representation for transcriptomic analyses in closely related Phoenix species. These results support the suitability of a reference-guided approach for comparative gene expression analysis between P. theophrasti and P. pusilla.

Table 1

Table 1: Summary of RNA-seq data.

Species-associated differential gene expression

To investigate transcriptional divergence between P. theo phrasti and P. pusilla, differential gene expression analysis was performed at the gene level. After filtering low-abundance tran scripts, a large number of genes exhibited significant differen t ial expression between the two species (adjusted p< 0.05). In total, 5,619 genes were up-regulated, and 5,901 genes were down-regulated in P. theophrasti relative to P. pusilla (Figure 1A), indicating extensive yet balanced transcriptional diver gence. Visualization of differential expression using a volcano plot revealed a broad distribution of fold changes accompa nied by strong statistical support (Figure 1B). More than 7,000 genes displayed large absolute log2fold changes together with extremely low adjusted p-values, suggesting that the observed expression differences are robust and not driven by marginal effects. Both up- and down-regulated gene sets spanned a wide range of functional annotations, indicating that multiple biolog ical processes contribute to species-specific expression profiles.

Figure 1: Differential gene expression between P. theophrasti and P. pusilla. (A) Bar plot showing the numbers of up- and down-regulated genes. (B) Volcano plot showing fold change vs significance.

Expression patterns of highly divergent genes

To visualize expression patterns among the most strongly differentiated genes, a heatmap was generated using variance stabilized expression values for the top differentially expressed genes ranked by fold change and statistical significance. Hierar chical clustering revealed a clear grouping of samples by spe cies, with replicates clustering closely together within each spe cies (Figure 2).

The distinct clustering patterns observed across species, cou pled with minimal variation within species, demonstrate that the transcriptional differences identified are consistent and re producible at the species level. These results further support the presence of stable and biologically meaningful gene expression divergence between P. theophrasti and P. pusilla leaf tissues.

Figure 2: Heatmap of top differentially expressed genes between P. theophrasti and P. pusilla.

Functional enrichment of differentially expressed genes

To assess the biological significance of the observed tran scriptional differences, functional enrichment analyses were conducted on the differentially expressed gene set. Gene Ontol ogy (GO) enrichment analysis revealed significant over-repre sentation of biological processes related to secondary metabol ic processes, ribosome structure and biogenesis, and responses to biotic and abiotic stimuli (Figure 3A). Enriched stress-related categories included responses to jasmonic acid, wounding, hy drogen peroxide, and antimicrobial compounds, highlighting differences in stress perception and signaling between species.

Notably, multiple ribosome-associated terms, such as struc tural constituent of ribosome and cytosolic ribosomal subunit, were among the most significantly enriched categories. This en richment suggests coordinated differences in translational ca pacity or regulation between P. theophrasti and P. pusilla leaves.

KEGG pathway enrichment analysis further corroborated these findings (Figure 3B). Significantly enriched pathways in cluded phenylpropanoid biosynthesis, flavonoid and diterpe noid biosynthesis, plant–pathogen interaction, MAPK signaling, and ribosome-related pathways. Together, these results indi cate that species-specific transcriptional divergence involves both core cellular processes and adaptive metabolic and signal ing pathways relevant to leaf physiology.

Figure 3: Functional enrichment of differentially expressed genes between P. theophrasti and P. pusilla. (A) Gene Ontology (GO) biological process enrichment of differentially expressed genes, showing significantly over-represented functional catego ries based on gene count and adjusted p-value. (B) Plant-relevant KEGG pathway enrichment of differentially expressed genes, high lighting pathways associated with metabolism, stress signaling, and defense responses.

Identification of trait-relevant candidate genes

To identify specific genes with potentially underlying func t ional divergence between species, differentially expressed genes were further examined at the individual gene level. High confidence candidate genes were selected based on annota t ion quality, expression magnitude, and statistical significance. These genes were grouped into trait-relevant functional catego ries, including abiotic stress signaling, oxidative stress and Re active Oxygen Species (ROS) regulation, defense and immunity, photosynthesis and light responses, cuticle and lipid metabo lism, and secondary metabolism (Table 2 & Figure 4).

Genes associated with stress signaling and kinase-mediated pathways were prominently represented, including receptor like kinases, calcium-associated signaling components, and heat shock proteins. Numerous genes involved in oxidative stress regulation, such as peroxidases, glutathione S-transferases, su peroxide dismutases, and oxidoreductases, exhibited strong dif ferential expression, suggesting species-specific modulation of redox homeostasis.

Defense- and immunity-related genes, including disease re sistance proteins, pathogenesis-related proteins, and NB-ARC domain-containing genes, also showed pronounced expres sion differences, consistent with enrichment of plant–pathogen interaction pathways. In addition, multiple photosynthesis related genes, including components of photosystem I and II, chlorophyll-binding proteins, and ATP synthase subunits, were differentially expressed, indicating divergence in photosynthet ic regulation at the transcript level.

Genes associated with cuticle formation, wax and lipid me tabolism, and secondary metabolic pathways, such as chalcone and terpene biosynthesis, were likewise differentially regulated between species. Collectively, these candidate genes provide mechanistic support for the enriched functional categories and suggest that transcriptional divergence between P. theophrasti and P. pusilla reflects coordinated differences in leaf metabo lism, stress responses, and physiological regulation.

Figure 4: Differentially Expressed Genes (DEGs) between P. theophrasti and P. pusilla. (A) Abiotic stress signaling. (B) Oxida t ive stress / ROS. (C) Defense and Immunity. (D) Photosynthesis / Light response. (E) Cuticle, wax, and lipid metabolism. (F) Second ary metabolism

Table 2

Table 2: Candidate genes associated with trait-relevant functional categories showing differential expression between P. theophrasti and P. pusilla.

Comparative transcriptome analysis of mature leaf tissues from P. theophrasti and P. pusilla revealed extensive and coor dinated transcriptional divergence between these two closely related palm species. Despite the use of a single reference genome (P. dactylifera), high alignment rates and consistent expression profiles across libraries indicate that the observed differences reflect genuine species-associated transcriptional variation rather than technical artifacts. Together, differential expression, functional enrichment, and trait-level gene analy ses point to divergence in multiple core biological processes, including stress signaling, metabolic regulation, photosynthe sis, and defense responses. The suitability of P. dactylifera as a reference genome for comparative analyses within the genus is supported by the availability of a high-quality genome assembly and extensive gene model annotation [1,2].

Species-level transcriptional divergence

The identification of more than 11,000 differentially ex pressed genes between P. theophrasti and P. pusilla highlights substantial transcriptional differentiation at the leaf level. The relatively balanced distribution of up- and down-regulated genes suggests widespread regulatory divergence rather than directional shifts confined to a single pathway. Similar mag nitudes of interspecific transcriptional divergence have been reported in comparative leaf transcriptome studies of closely related plant species, even when overall genome structure re mains highly conserved [7,34].

Clear species-specific clustering in the expression heatmap further supports the reproducibility and coherence of these transcriptional differences. Comparable patterns of clustering by species have been observed in RNA-seq studies of wild and cultivated monocots, including rice, maize, and other grasses, where leaf transcriptomes reflect both developmental pro gramming and environmental adaptation [19,37]. In palms, transcriptome-level differentiation has been less extensively explored, making the present study an important contribution to understanding molecular divergence within the genus Phoe nix [46]. Similar leaf-level transcriptome coherence has also been reported in date palm and royal palm under salinity and drought stress conditions [14,18,35,43].

Divergence in translational capacity and cellular processes

One of the findings of the functional enrichment analysis was the significant over-representation of ribosome-related GO terms and KEGG pathways. Enrichment of genes encoding ribo somal proteins and translational machinery suggests differenc es in protein synthesis capacity or regulation between the two species. Differential regulation of ribosomal genes has been reported in multiple plant comparative transcriptomics studies and is often associated with variation in growth rates, metabolic activity, or stress responsiveness [13,14,22,44].

In leaves, translational regulation plays a central role in co ordinating photosynthesis, stress responses, and metabolic f luxes. The observed divergence in ribosome-associated genes therefore likely reflects broader differences in cellular regula t ion and resource allocation between P. theophrasti and P. pusilla, rather than isolated effects on individual pathways. Enrichment of ribosome-related transcripts has similarly been observed in reference-guided transcriptome analyses of date palm tissues [31].

Divergence in stress signaling and redox regulation

Genes associated with abiotic stress signaling, oxidative stress, and Reactive Oxygen Species (ROS) regulation constitut ed a major component of the differentially expressed gene set. A number of receptor-like kinases, calcium-associated signaling proteins, and heat shock proteins exhibited strong species-spe cific expression patterns. Such signaling components are cen tral to plant perception and integration of environmental cues, particularly in leaf tissues exposed to fluctuating light, tempera ture, and water availability [29,47].

Similarly, the prominent representation of peroxidases, glu tathione S-transferases, superoxide dismutases, and oxidore ductases indicates divergence in redox homeostasis between species. Comparative studies in monocots have shown that variation in ROS-related gene expression often underlies dif ferences in stress tolerance and acclimation capacity [15]. Date palm transcriptome studies under heat, drought, salinity, and cadmium stress have consistently reported differential regula t ion of antioxidant enzymes, stress-responsive transcription factors, and signaling components in leaf tissues [18,33].

Defense- and immunity-related transcriptional differences

Defense-related genes, including disease resistance proteins, pathogenesis-related proteins, and NB-ARC domain-containing genes, were among the most strongly differentiated between P. theophrasti and P. pusilla. Enrichment of plant–pathogen inter action pathways further support divergence in immune-related transcriptional programs. Similar patterns have been observed in comparative transcriptome analyses across plant species, where immune genes are often among the most rapidly evolv ing and differentially regulated components of the genome [17,20].

Leaf tissues represent a primary interface between plants and their biotic environment, and variation in immune gene expression may reflect differences in pathogen pressure or eco logical context experienced by the two species. While function al assays would be required to confirm adaptive significance, the transcriptional patterns observed here suggest that defense signaling constitutes a key axis of divergence within Phoenix. Genome-wide analyses of date palm have also reported expan sion and diversification of resistance gene families, underscor ing the importance of immune regulation in palms [2].

Photosynthesis and leaf metabolic regulation

Differential expression of genes encoding photosystem com ponents, chlorophyll-binding proteins, ATP synthase subunits, and photosystem assembly factors indicate divergence in pho tosynthetic regulation between species. Comparative transcrip tomics in leaves has repeatedly demonstrated that even subtle differences in photosynthesis-related gene expression can in f luence photosynthetic efficiency, light utilization, and overall carbon metabolism [11,25].

In addition to photosynthesis, genes associated with cuticle formation, wax biosynthesis, lipid metabolism, and secondary metabolism were differentially expressed. Cuticle and wax com ponents play critical roles in water retention, protection from environmental stress, and leaf surface interactions [36]. Diver gence in these pathways has been reported in monocots occu pying contrasting environments and is often linked to ecological specialization [32]. Leaf transcriptome analyses in date palm have similarly shown coordinated regulation of photosynthesis- and cuticle-related genes under salinity stress [18].

Secondary metabolism and species-specific chemical pro f iles

The enrichment of phenylpropanoid, flavonoid, and terpene biosynthesis pathways, together with differential expression of chalcone synthases, terpene synthases, and related enzymes, suggests divergence in secondary metabolic capacity between P. theophrasti and P. pusilla. Secondary metabolites contribute to defense, stress tolerance, and interactions with herbivores and microbes, and their regulation is frequently species-specific [10,38]. Comparative transcriptomic studies in monocots and woody plants have shown that variation in secondary metabo lism often reflects both evolutionary history and environmental adaptation [12,26]. Enrichment of secondary metabolic path ways has also been reported in date palm transcriptomes under abiotic stress, suggesting conserved regulation of these path ways across palms [12].

Methodological considerations

Because RNA samples were pooled for each species, the present study captures species-level transcriptional differenc es rather than within-species variation. Similar pooled sample strategies have been successfully employed in comparative tran scriptomics when the primary objective is interspecific compari son [8,30]. The use of P. dactylifera as a reference genome rep resents a pragmatic and widely accepted strategy in non-model plant transcriptomics. The high alignment rates observed here support the validity of this approach, although species-specific transcripts absent from the reference genome may not be fully captured. Comparable reference-guided approaches have been successfully applied in multiple date palm transcriptome stud ies across tissues and stress conditions [18,31].

This comparative transcriptome analysis reveals extensive and coordinated transcriptional divergence between P. theo phrasti and P. pusilla palms. Differences in translational ma chinery, stress signaling, redox regulation, defense responses, photosynthesis, and secondary metabolism collectively suggest that species-specific transcriptional programs underlie func t ional differentiation within the genus Phoenix. By integrating differential expression, functional enrichment, and trait-focused gene analysis, this study provides a comprehensive framework for understanding transcriptome divergence in palms and es tablishes a foundation for future ecological, physiological, pop ulational and evolutionary investigations.

Author contributions

VMG and MNR analyzed, reviewed, and drafted the manu script. MNR conceptualized the experiments, edited the manu script, and contributed to the logistical support for the execu t ion of the experiment.

Conflict of interest

Authors in this research work declare no conflict of interest.

Acknowledgements

This research was supported through appropriated funds from the United States Department of Agriculture - Agricul tural Research Service (USDA-ARS) under Project Number 603813210-004-000-D. Various biological technicians who maintain and help with sampling palms are highly acknowl edged.