*4.4. RNA-Seq Data Analysis*

After removing adapter sequences, ploy-N-containing reads, and low-quality sequences from the raw data, the remaining clean reads were mapped to the reference genome via Hisat2 (version 2.1.0) software [53]. The expression levels of each gene were calculated and normalized as fragments per kilobase million (FPKM) measurements. Functional annotations of the identified transcripts were then aligned against various databases, including NR [54], Swiss-Prot [55], protein family (Pfam) [56], Gene Ontology (GO) [57], and Kyoto Encyclopedia of Genes and Genomes (KEGG) [58]. Differential gene expression analysis across varieties was performed using the DESeq2 package (version 1.16.1) [59], and genes with an adjusted *p*-value < 0.05 and absolute log2 fold changes (FC) ≥ 1 were set as the screening criteria for significantly differential expression. Based on the identified DEGs, three types of GO categories were obtained, including biological processes (BP), cellular components (CC), and molecular function (MF). Additionally, the functional enrichment analysis of GO and KEGG was implemented to identify which genes were significantly enriched. GO terms and KEGG pathway with a threshold of false discovery rate (FDR) < 0.05 were considered as the significance level.
