The genomic and annotation data utilized in this study were obtained from two sources. First, six Rhus genomes assembled by our laboratory were included, achieving gold-standard reference genome quality to ensure the reliability and completeness of this research. Second, publicly available genomes of three additional Anacardiaceae species Mangifera indica, NCBI accession GCA_011075055.1 (Wang et al., 2020c), Pistacia vera, GCA_008641045.1 (Zeng et al., 2019), and Anacardium occidentale (Phytozome13, https://phytozome-next.jgi.doe.gov) were incorporated to broaden the phylogenetic scope (Goodstein et al., 2011). To provide an evolutionary framework for comparative analyses, the Arabidopsis thaliana genome (TAIR version 10) was included as an outgroup (Lamesch et al., 2011). This integrated dataset enabled a comprehensive investigation of genomic diversity and evolutionary dynamics within Anacardiaceae ( Supplementary Table S1 ).

Genome collinearity was investigated through JCVI (v1.2.7) using protein-based reciprocal best-hit alignments to detect conserved syntenic blocks (Tang et al., 2008). Putative orthologous relationships were established through pairwise genome comparisons, with aligned regions representing ancestral genomic linkages. Whole-genome duplication (WGD) chronology and speciation divergence times were resolved via WGDI (v0.6.1) analysis (Sun et al., 2022), incorporating calculation of nonsynonymous (Ka) and synonymous (Ks) substitution rates through the NG86 method (Nei and Gojobori, 1986). This integrative approach enabled temporal reconstruction of lineage-specific WGD events and interspecies divergence patterns.

Full-length long terminal repeat retrotransposons (LTR-RTs) were identified through complementary de novo prediction using LTRharvest (v1.5.10) (Ellinghaus et al., 2008) and LTR_FINDER_parallel (v1.1) (Ou and Jiang, 2019). Consensus predictions were integrated and refined through LTR_retriever (v2.9.0) (Ou and Jiang, 2018), with LTR insertion timing calculated using a lineage-specific substitution rate of 7×10^-9 (Ossowski et al., 2010).

Genome-wide SNP identification was performed through MUMmer (v4.0.0beta2) alignment using stringent filtering parameters (-r -q -l 1000) to ensure reciprocal best matches (Marçais et al., 2018). The show-snps utility was implemented with -ClrTH thresholds to detect single nucleotide variants and small indels (< 100 bp). Functional impacts of identified polymorphisms were annotated using SnpEff (v5.0e) through comprehensive variant effect prediction (Cingolani et al., 2012). Visualization of the interchromosomal distribution of structural variations in Rhus species using SyRI (Goel et al., 2019).

To investigate the evolutionary relationships within Anacardiaceae, we performed phylogenetic reconstruction incorporating six Rhus species from this study along with three additional Anacardiaceae taxa (P. vera, M. indica, A. occidentale) and the outgroup species A. thaliana. Single-copy ortholog identification was achieved through OrthoFinder (v2.5.4) analysis (Emms and Kelly, 2019), employing the longest transcript isoforms from ten plant species as input data. Multiple sequence alignment of conserved orthologous sequences was executed using MUSCLE (v5.1.linux64) (Edgar, 2022), followed by identification of phylogenetically informative regions through Gblocks (v0.91b) for site selection (Castresana, 2000).

Optimal substitution model selection via Prottest (v3.4.2) identified Blosum62+I+G+F as the most appropriate amino acid substitution matrix (Darriba et al., 2011). Maximum likelihood phylogenetic inference was subsequently conducted using RAxML (v8.2.12) with 200 bootstrap replicates to assess nodal support (Stamatakis, 2014). Temporal calibration of divergence events was implemented in MCMCTREE (PAML v4.9j) (Yang, 2007), incorporating a key fossil-derived calibration point (Rhus-P. vera split at ~48 mya) from TimeTree and prior paleobotanical studies (Moran, 1989; Kumar et al., 2017).

Gene family dynamics were analyzed through CAFE (v5.0) (Mendes et al., 2021), leveraging OrthoFinder-derived gene clusters to detect significant expansions/contractions (conditional P < 0.05 threshold). Functional enrichment analysis of dynamically evolving gene families was conducted using the R package clusterProfiler by performing Gene Ontology (GO) term and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses (Ogata et al., 1999; Ashburner et al., 2000; Yu, 2024). KEGG pathway annotations were retrieved from the official database (http://www.kegg.jp/kegg/pathway.html).

To systematically identify and characterize defense-associated genes, we implemented an integrated annotation framework combining structural domain verification and functional ontology analysis. The WRKY transcription factors were rigorously identified by detecting the presence of the PF03106 WRKY domain through Pfam profiles. Specifically, the Pfam database integrated within the InterPro platform (https://www.ebi.ac.uk/interpro/) was queried using Hidden Markov Model (HMM) profiles with a stringent E-value cutoff of 1e-5 to ensure high-confidence annotation of these conserved protein domains (Blum et al., 2024). NLR-type resistance genes were classified based on the presence of characteristic nucleotide-binding domains (PF00931/NB-ARC, PF01582/TIR, PF05659/RPW8) coupled with leucine-rich repeat motifs (PF00560/LRR_1, PF07725/LRR_8, PF13306/LRR_9, PF13855/LRR_6). Functional annotation was extended through GO enrichment analysis targeting seven phytohormone signaling pathways: auxin (GO:0009734), cytokinin (GO:0009736), abscisic acid (GO:0009738), brassinosteroid (GO:0009742), jasmonic acid (GO:0009867), ethylene-activated signaling pathway (GO:0009873), and strigolactone (GO:0010476) responses. Additionally, selective pressure analysis was performed using PAML to calculate the dN/dS ratios of defense-related genes relative to background genes.

To investigate the evolutionary relationships of defense-associated genes in Anacardiaceae, multiple sequence alignment of WRKY proteins identified in Rhus and Arabidopsis was performed using MUSCLE (v5.1.linux64). A neighbor-joining phylogenetic tree was subsequently constructed with 1,000 bootstrap replicates in MEGA11 to assess node reliability (Koichiro et al., 2021). Conserved motifs were predicted using the MEME Suite (v5.5.7) with default parameters (Bailey et al., 2015), while functional domains were annotated via the NCBI Conserved Domain Database (CDD) (Marchler-Bauer et al., 2014). Gene structure features (coding sequences and untranslated regions) and chromosomal localization patterns were visualized using TBtools (v2.225) to integrate genomic architecture and positional context (Chen et al., 2023).

We retrieved 2-kb upstream promoter regions of R. chinensis WRKY family genes using genome annotation data. Putative cis-regulatory elements were systematically characterized through computational analysis using PlantCARE (https://bioinformatics.psb.ugent.be/webtools/plantcare/html/), with particular emphasis on stress-responsive motifs. Core promoter elements, such as the TATA-box and CAAT-box, were excluded from subsequent analyses to specifically focus on defense-related regulatory sequences. Visualization of spatial element distribution was performed using TBtools with customized annotation parameters (Chen et al., 2023). For protein interaction network inference, a homology-based approach was implemented. All annotated R. chinensis protein sequences were compiled into a local BLAST database. WRKY family proteins served as query sequences for iterative BLAST searches with E-value < 1e-5 and sequence identity >= 30%. The network topology was constructed using Gephi (Bastian et al., 2009), in which nodes represent individual WRKY proteins and edges denote significant sequence homology, indicating potential functional associations.

Following specimen acquisition, tissues were rapidly preserved through cryogenic freezing to maintain biomolecular stability. Nucleic acid isolation was performed employing commercial plant RNA extraction kits (Tiangen Biotech) with subsequent enzymatic purification using recombinant DNase I (Takara Bio) to eliminate genomic DNA interference. RNA integrity verification was conducted through nucleic acid electrophoresis utilizing 1% agarose matrices impregnated with fluorescent nucleic acid stain, supplemented by dual-platform validation of nucleic acid purity and concentration through spectrophotometric (NanoDrop 2000C) and microfluidic (Agilent Bioanalyzer 2100) quantification systems. For Illumina sequencing, poly-A enriched transcripts were converted into amplifiable cDNA libraries through reverse transcription employing ultra-high-fidelity cDNA synthesis reagents (NEBNext series, NEB). The resulting PCR products were purified using the AMPure XP system, and the library quality was assessed on the Agilent Bioanalyzer 2100 system.

To ensure data reliability, three biological replicates were included for each sample. Raw sequencing data underwent comprehensive quality evaluation through FastQC (v0.12.1) prior to pre-processing (Andrews, 2010). Adapter trimming and read quality enhancement were implemented via Trim Galore (v0.6.7) (Krueger et al., 2021), followed by alignment of processed reads to the reference genome using HISAT2 (v2.2.1) with standard parameters (Kim et al., 2015). Transcript assembly was performed through StringTie (v2.1.0) (Pertea et al., 2015), with subsequent quantification of gene expression levels expressed as fragments per kilobase of transcript per million fragments mapped (FPKM). Differential expression analysis was executed through DESeq2 (v1.38.3) (Love et al., 2014), employing a stringent significance threshold (P < 0.01) coupled with a |log2 fold change| ≥ 1 for candidate gene identification.

The statistical details of analysis applied in this paper are provided alongside in the results and methods section. Statistical analyses were performed in R 4.0.0 (R Core Team, 2014).

The Anacardiaceae family, which includes economically and ecologically important genera such as Rhus, Mangifera, and Pistacia, exhibits substantial genomic diversity and evolutionary complexity. To explore these variations, we used R. punjabensis as the reference genome for SNP calling across five other Rhus species (R. chinensis, R. chinensis var. roxburghii, R. typhina, R. potaninii, and R. glabra). Analysis of the chromosomal distribution revealed distinct compartmentalization, with the B compartment concentrated in centromeric regions, characterized by higher transposon density, while the A compartment was primarily distributed along chromosome arms, enriched with gene-dense regions ( Figure 1B ). This spatial organization reflects a high-quality genome assembly, reinforcing the reliability of SNP identification and interspecies comparisons. To further expand our comparative framework, we included three economically important species from the Anacardiaceae family: P. vera, M. indica, and A. occidentale, providing a broader perspective on genomic evolution within the family. As shown in Figure 1C , genome sizes varied significantly, ranging from 324.4 Mb in R. potaninii to 671.1 Mb in P. vera, with the number of genes spanning 20.6 k to 41.3 k ( Supplementary Table S2 ). Notably, when P. vera was excluded, genome sizes of the remaining species showed relatively minor differences, indicating that large-scale expansions may be lineage-specific rather than characteristic of the entire family. Synteny analysis between R. chinensis and M. indica ( Figure 1D ) revealed a clear 1:2 syntenic relationship, suggesting that M. indica underwent a WGD event followed by chromosomal rearrangements, resulting in its 20-chromosome configuration. In contrast, Rhus species displayed a conserved chromosomal structure without evidence of recent WGD, implying that their diversification was primarily driven by smaller-scale duplications, SNP accumulation, and transposable element activities. These findings highlight not only the high assembly quality of the Rhus genomes but also the distinct evolutionary trajectories among major Anacardiaceae lineages.

Phylogenetic analysis based on whole-genome sequencing and orthologous gene families from nine representative Anacardiaceae species, with A. thaliana as an outgroup, enabled the reconstruction of evolutionary trajectories and divergence times. The resulting phylogenetic tree revealed that Rhus species form a distinct clade closely related to Pistacia. Molecular dating suggests that Anacardiaceae originated approximately 97 million years ago (Mya), with Rhus diverging from its common ancestor around 48 Mya ( Figure 2A ). This split likely reflects geographical isolation and subsequent ecological adaptation of Rhus across temperate and subtropical regions. Notably, significant gene family expansions were observed in P. vera, M. indica, and A. occidentale, potentially driven by WGD events and transposable element proliferation. Prior studies indicate that M. indica and A. occidentale share a common WGD event (Mango Genome Consortium et al., 2021), whereas synteny analysis between R. chinensis and M. indica demonstrated that Rhus did not undergo an independent WGD ( Figure 2B ). Comparative genomic analysis revealed that Rhus shares 13,369 gene families with three other Anacardiaceae species while possessing 1,652 unique gene families ( Figure 2C ). Further syntenic analysis with the ancestral eudicot karyotype (AEK) confirmed that Rhus experienced only the gamma duplication event shared among eudicots, without any lineage-specific WGD ( Figure 2D ). In total, 23,530 genes from R. chinensis were clustered into 16,411 gene families, of which 81.5% were conserved across Rhus and three other Anacardiaceae species, while 10.0% represented Rhus-specific families ( Figures 2C, E ). These unique gene families may contribute to the species-specific traits and ecological adaptability of Rhus.

We conducted a comprehensive comparative analysis of structural variations (SVs) across Rhus species using R. chinensis as the reference genome. The analysis identified substantial differences in the number and types of genetic variations, including single nucleotide polymorphisms (SNPs), insertions/deletions (InDels), and SVs such as presence-absence variations (PAVs), inversions, and translocations. Among these species, R. chinensis var. roxburghii exhibited the highest number of SVs, with a total of 6,099 events, followed closely by R. glabra (5,355) and R. typhina (5,218). In terms of specific variation types, R. potaninii and R. glabra displayed the greatest number of deletions, with 143,237 and 143,254 events, respectively. Notably, the number of insertions was also prominent in R. potaninii (132,100), while R. punjabensis demonstrated the highest count of translocations (1,131). Comparative analysis illustrated clear lineage-specific patterns of SV distribution, as reflected in the contrasting SV profiles between R. chinensis var. roxburghii and R. chinensis ( Supplementary Figure S1 ), R. punjabensis and R. potaninii ( Supplementary Figure S2 ), and R. typhina and R. glabra ( Supplementary Figure S3 ). These variations may underlie phenotypic divergence, ecological adaptation, and evolutionary differentiation among Anacardiaceae species. The identification of such structural variations provides valuable insights into the genomic mechanisms contributing to species-specific traits and adaptive capacities within Anacardiaceae.

Furthermore, transposable elements (TEs) have played a pivotal role in shaping genome architecture within Anacardiaceae. In contrast to M. indica, P. vera does not exhibit lineage-specific genome duplications, although it shares the ancient gamma duplication event common to eudicots. Instead, its larger genome size is primarily driven by substantial TE expansion, with TEs occupying a remarkable 70.7% of its genome and LTR-RTs accounting for 46.75%, rather than typical WGD observed in other species (Zeng et al., 2019). The insertion time distribution analysis of LTR-RTs in various Anacardiaceae species reveals multiple waves of transposition, potentially triggered by environmental shifts that facilitated adaptive evolution ( Figure 2F ). In Rhus, the insertion time analysis of intact LTR-RTs demonstrated that R. chinensis and R. chinensis var. roxburghii exhibit similar insertion patterns ( Supplementary Figure S4 ), with recent transposition bursts observed in both lineages. Likewise, R. typhina and R. glabra showed comparable insertion dynamics ( Supplementary Figure S5 ). In contrast, R. punjabensis and R. potaninii displayed distinct LTR insertion timelines, with R. potaninii experiencing a transposition burst approximately two million years ago, while R. punjabensis showed more recent LTR insertions around one million years ago ( Supplementary Figure S6 ). The chromosomal distribution of retrotransposons suggests a potential role in enhancing genomic plasticity and adaptive capacity.

Gene family expansion, a critical driver of adaptive evolution (Guo, 2013), was analyzed using the CAFE program across nine Anacardiaceae species. Notably, P. vera, M. indica, and A. occidentale exhibited significant expansions of 2,609, 2,048, and 3,927 gene families, respectively, with expansions notably outpacing contractions, suggesting that these species undergo adaptive evolution primarily through the expansion of specific gene families. Functional enrichment analyses revealed defense-related pathways as central to these expansions. For P. vera ( Figure 3 ; Supplementary Tables S3 , S4 ), GO terms such as “defense response to bacterium,” “xenobiotic detoxification by transmembrane export across the plasma membrane,” and “protein phosphorylation” were prominently enriched ( Figure 3A ). Correspondingly, KEGG pathways included “plant−pathogen interaction” and “biosynthesis of various plant secondary metabolites” ( Figure 3B ), highlighting their involvement in redox regulation, phytohormone-mediated defense signaling, and defense-associated metabolic processes. In Anacardiaceae, enriched GO terms included “defense response to bacterium,” “defense response,” and “signal transduction,” reflecting defense-related biological processes ( Supplementary Figure S7 ; Supplementary Tables S5 , S6 ). Similarly, M. indica showed significant enrichment in “regulation of defense response to fungus” and “response to auxin,” suggesting that the expanded gene families are predominantly associated with defense-related functions ( Supplementary Figure S8 ; Supplementary Tables S7 , S8 ). For A. occidentale, KEGG pathways such as “Plant hormone signal transduction” were highlighted, indicating its involvement in hormone-mediated defense responses ( Supplementary Figure S9 ; Supplementary Tables S9 , S10 ).

WRKY transcription factors (TFs) are pivotal regulators of plant defense mechanisms, orchestrating responses to biotic and abiotic stresses (Jiang et al., 2017). In R. chinensis, a comprehensive analysis was performed to investigate the phylogenetic relationships, conserved domains, and structural organization of WRKY TFs. Based on the conserved WRKY domain PF03106, a total of 55 WRKY genes were identified in R. chinensis. To elucidate their evolutionary relationships, 15 representative WRKY sequences from A. thaliana were retrieved from the TAIR database and aligned with the 55 R. chinensis WRKY sequences for phylogenetic analysis ( Supplementary Figure S10 ; Supplementary Table 11 ). Structural analysis of the WRKY proteins revealed the presence of 10 conserved motifs distributed across the 55 members. These motifs exhibited distinct patterns among different WRKY subgroups, suggesting divergent evolutionary trajectories and specialized physiological roles ( Figure 4 ; Supplementary Table S12 ). All R. chinensis WRKY proteins contained at least one WRKY superfamily domain. Additionally, several members harbored other conserved domains, such as plant-specific zinc finger motifs and bZIP superfamily elements, which may enhance transcriptional regulation under pathogen invasion and environmental challenges ( Figure 4 ). Gene structure analysis further highlighted substantial variation in intron-exon arrangements among WRKY members. Some genes exhibited multiple introns, indicative of potential alternative splicing events that could contribute to functional diversity and regulatory complexity.

The functional diversity of genes is closely associated with the composition and abundance of cis-regulatory elements in their promoter regions. To investigate the potential regulatory mechanisms of WRKY transcription factors, we systematically analysed 2 kb upstream sequences of WRKY genes, with particular emphasis on cis-elements associated with abiotic stress adaptation, biotic stress responses, and defense metabolism regulation ( Figure 5A ). Our analysis revealed three predominant regulatory elements: MeJA-responsiveness (involved in jasmonate signaling), abscisic acid (ABA) responsiveness (central to drought and osmotic stress adaptation), and salicylic acid (SA) responsiveness (critical for systemic acquired resistance against pathogens), which emerged as the most prominently represented stress-related motifs ( Figure 5B ). These elements coordinate plant responses to environmental challenges through phytohormone-mediated signalling pathways. Furthermore, the promoter regions contained multiple specialized cis-regulatory modules, including elements governing low-temperature adaptation (e.g., dehydration, low-temperature, and salt stress-responsive motifs), hypoxia-inducible enhancer-like sequences, and MYB-binding sites regulating both drought-inducible gene expression and flavonoid biosynthesis. Notably, wound-responsive elements and pathogen defense motifs (e.g., elicitor-mediated activation elements) were selectively retained in specific WRKY promoters, implying specialized adaptation for localized stress recognition. The homology-based protein interaction network, constructed using sequence similarity thresholds, revealed a striking pattern of intra-family connectivity among WRKY proteins ( Figure 5C ). Nodes representing WRKY members exhibited dense clustering, with extensive edge formation signifying robust sequence homology-driven associations. Notably, no significant interactions were detected between WRKY proteins and non-WRKY proteins, suggesting functional modularity or selective evolutionary constraints that prioritize intramolecular cooperation within the WRKY family.

Chromosomal mapping revealed that WRKY, NLR, and phytohormone-related genes are unevenly distributed across the chromosomes of R. chinensis. Notably, these genes exhibit clustered arrangements in specific chromosomal regions ( Figure 6A ; Supplementary Figure S11 ). Among them, NLR genes are particularly enriched on chromosomes 4 and 12 ( Supplementary Figure S12 ), suggesting that these loci may play crucial roles in adaptive defense mechanisms. In R. chinensis, genome-wide identification and characterization of these defense-associated genes indicate their extensive presence across multiple chromosomes, with pronounced clustering at certain loci. This non-random distribution pattern implies potential hotspots for evolutionary adaptation, possibly driven by selective pressures to enhance pathogen defense. To further understand the evolutionary dynamics of these genes, we performed evolutionary analysis using the dN/dS ratio (nonsynonymous to synonymous substitution ratio). The analysis identified evidence of positive selection in these defense-related gene families (Student’s t-test, P < 0.05; Figure 6B ), indicating that mutations in these loci are likely beneficial for improving defense responses against biotic stresses. These findings underscore the evolutionary importance of defense-associated genes in shaping the adaptive landscape of Rhus species.

In this study, the synteny relationships of defense-related genes were analyzed among six Rhus species and three additional members of the Anacardiaceae family, namely M. indica, P. vera, and A. occidentale. High-resolution synteny maps revealed both the conservation and specificity of defense-related gene loci, with particular emphasis on regions associated with the WRKY gene family. The analysis primarily focused on syntenic clusters of WRKY-related genes, as illustrated in Figure 6C . Notably, conserved WRKY gene fragments were identified at the terminal or proximal regions of chromosomes 2, 3, and 7 across the six Rhus species, indicating strong syntenic relationships. Typically, these segments corresponded to a single copy in P. vera, whereas M. indica and A. occidentale maintained two copies, suggesting lineage-specific duplication events. Interestingly, certain deviations were observed. For example, on chromosome 7, the two North American Rhus species (R. typhina and R. glabra) exhibited specific loss of syntenic segments compared to the four East Asian Rhus species. Conversely, these syntenic regions were well-preserved in the more distantly related P. vera, indicating potential species-specific adaptive evolution.

To validate the functional relevance of defense-associated genes identified in Rhus species, we investigated their transcriptional responses to infestation by gall-forming aphids, obligate parasites that induce tannin-rich galls (Wei et al., 2022). Samples were collected from four tissue-time point combinations: leaf-wing before infection (WBI), leaf-wing after 7 days of infection (WAI-7d), leaf before infection (LBI), and leaf after 7 days of infection (LAI-7d). RNA-seq analysis revealed that biotic stress significantly upregulated 31 WRKY transcription factor genes in leaves (LAI-7d vs. LBI; fold change >2; Figure 6D ; Supplementary Tables S13 , S14 ). Moreover, WRKY genes in leaf-wings (WAI-7d vs. WBI) exhibited a broad upregulation trend across the gene family, with expression levels increasing in the majority of members, though these changes did not reach statistical significance. Similarly, defense-related genes, including NLRs and phytohormone signaling components, exhibited pronounced expression shifts post-stress ( Supplementary Figure S13 ). These results demonstrate that defense-related gene families under positive selection are functionally responsive to biotic stress.