The basic information on genome assembly and genes was obtained from the genome and annotation files. Following the published genome annotation pipeline (Haas et al., 2008, Holt and Yandell, 2011, Thibaud et al., 2016, Hoff et al., 2019), structural annotation was performed on 15 genomes by June 2024. A species-specific repeat library was built using RepeatModeler 2.03 (https://www.repeatmasker.org/) and combined with the RepBase 21.12 database (Bao et al., 2015) to comprehensively identify repetitive elements across the genome using RepeatMasker 4.13 (Flynn et al., 2020) with default parameters. Four tools were used for de novo gene prediction: GeneMark-ET 4.38 (Borodovsky and McIninch, 1993) with default parameters; SNAP 7.0 (Korf, 2004) with default parameters; AUGUSTUS 3.4.0 (Stanke et al., 2006) --species --gff3=on; and BRAKER 2.1.6 (Bruna et al., 2021) run via braker.pl --genome genome.fa --bam rnaseq.bam --species=eucalyptus. Homology-based predictions were performed with GenomeThreader 1.7.3 (Myburg et al., 2014) with default parameters by aligning protein sequences from Arabidopsis thaliana, Oryza sativa, Populus trichocarpa, and Eucalyptus (Myburg et al., 2014). For transcriptome-based evidence, StringTie 2.2.1 (Pertea et al., 2015) was utilized for the de novo assembly of RNA-Seq data from Eucalyptus. Trinity 2.8.5 (Haas et al., 2003) with --genome_guided_bam rnaseq.bam was employed for genome-guided transcriptome assembly. TransDecoder 5.5.0 (Haas et al., 2013) was then applied to predict the coding sequences from the assembled transcripts. RNA-Seq data were integrated using PASA 2.5.2 (Haas et al., 2008) via Launch_PASA_pipeline.pl -c alignAssembly.config -R -g genome.fa -t transcripts.fa --ALIGNERS blat,gmap to construct a comprehensive transcriptome database. All evidence types—de novo predictions, homology-based annotations, and transcripts—were integrated in MAKER3 3.1.3 (Cantarel et al., 2008) using maker -base eucalyptus_annotation -genome genome.fa -est transcripts.fa -protein proteins.fa -rm_gff repeatmasker_out.gff3 -pred_gff predictions.gff3 to optimize gene model prediction. EVM 2.1.0 (Haas et al., 2008) --weights weights.txt --gene_predictions denovo_preds.gff3 --transcript_alignments transcript_out.gtf --protein_alignments proteins.gff3 was subsequently used to generate weighted consensus gene models. A final round of annotation updates was performed using PASA 2.5.2 (Haas et al., 2008) with -c alignAssembly.config -A -g genome.fasta -t all_transcripts.fasta –annots EVM.gff3 to refine gene structures and predict alternative splicing events and 5’ and 3’ untranslated regions. Gene function prediction was conducted for all 46 genomes. Gene functional annotation was performed by querying multiple databases, including Swiss-Prot (Apweiler et al., 2004), NCBI-NR (Pruitt et al., 2007), Pfam (Mistry et al., 2021), Gene Ontology (GO) (Ashburner et al., 2000), and Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2004).

  Non-coding RNA prediction was conducted for all 46 genomes. Candidate lncRNAs were identified based on full-length transcript assemblies generated with the Iso-Seq data. Transcripts with a length ≥200 nucleotides and containing at least 2 exons were initially chosen as lncRNAs. These were further filtered using computational tools with coding potential assessment capabilities, including CNCI 2.0 (Sun et al., 2013), CPAT 3.0.0, CPC2 0.1 (Kang et al., 2017), PLEK 2 (Li et al., 2024a), and Pfam 35.0 (Mistry et al., 2021), with default parameters, to reliably distinguish non-coding transcripts from protein-coding transcripts. tRNAs were predicted using tRNAscan-SE 2.0.12 (Lowe and Eddy, 1997) with default parameters. snRNAs and miRNAs were identified by aligning sequences against the Rfam database (Kalvari et al., 2018) through Infernal 1.1.4 (Nawrocki and Eddy, 2013) with cmscan --cut_ga --rfam --nohmmonly --tblout infernal.tblout Rfam.cm genome.fa. rRNAs were annotated using BLAST 2.14.0 (Kent, 2002) by aligning the eucalyptus genome against known rRNA sequences to accurately locate rRNA regions.

  Orthologous and paralogous genes were identified in 39 species (45 genomes) using OrthoFinder 2.5.4 (Emms and Kelly, 2019) with parameters -M msa -S diamond -T fasttree. Collinear blocks were identified using MCScanX 1.1.11 (Wang et al., 2012) with -e 1e-5 -s 5 -m 25. A phylogenetic tree was constructed using MEGA-X 10.2.6 (Kumar et al., 2018) applying the maximum likelihood method and 1000 bootstrap replicates based on single-copy orthologous gene families among eucalyptus species.

  After quality control using fastp 0.20.0 (Chen et al., 2018) with default parameters, the reads were aligned to the reference genome using HISAT2 2.0.5 (Kim et al., 2019) with --dta –new-summary --rna-strandness RF. Subsequently, Samtools 1.13 (Li et al., 2009) was further used to convert SAM format into BAM format. Gene expression was quantified using StringTie 1.3.6 (Pertea et al., 2015) with -e -B -G ref_annotation.gtf. For miRNA-seq analysis, sRNAminer 1.1.2 (Li et al., 2024b) with default parameters was employed in miRNA identification, target prediction, and expression analysis. A co-expression network was constructed using WGCNA 1.70.3 (Langfelder and Horvath, 2008). The TPM expression matrix was log2-transformed and quantile-normalized across samples using the preprocessCore package. Genes with zero variance or expressed in fewer than 50% of samples were excluded from the analysis. The network was built with the parameters power = 12, minModuleSize = 30, and mergeCutHeight = 0.25 to measure the correlation of gene expression patterns across samples using a weighted correlation matrix and to group genes into modules based on their co-expression patterns. The resulting co-expression network was visualized using Cytoscape 3.6.1 (Su et al., 2014) with default layout settings.

  For ChIP-seq, DAP-seq, and DNase-seq, Trimmomatic 0.36 (Bolger et al., 2014) with default parameters was used to remove adapter sequences and filter out low-quality reads. The clean data obtained from germplasm resources were aligned to the reference genome assemblies using Bowtie2 2.3.2 (Langdon, 2015) with default parameters. Picard 2.19 with REMOVE_DUPLICATES=true was used to remove PCR-duplicated reads, and peaks were called using the callpeak module of MACS2 2.1.2 with -nomodel -q 0.05 --extsize 200 --keep-dup all -B --call-summit (Feng et al., 2012).

  After using Trimmomatic 0.36 (Bolger et al., 2014) with default parameters to control the raw data quality, BWA 0.7.12 (Jo and Koh, 2015) was used to map the clean data to the reference genome. Samtools 1.13 (Li et al., 2009) was then used to transform and sort the mapping results, and GATK 4.1.5.0 (Summa et al., 2017) with parameters QUAL < 30.0, MQ < 50.0, QD < 2 was used to implement a unified genotype approach for SNP and InDel calling.