Genome-Wide Identification and Expression Analysis of Glutathione S-Transferases (GSTs) in Crocus sativus (Saffron) Reveals Their Roles in Hormone-Mediated Stress Responses

1. Introduction

Glutathione S-transferases (GSTs; EC 2.5.1.18) are phase II metabolic isozymes common in fungi, mammals, and plants [1]. The first report of GSTs in plants dates back to 1970, when research revealed that GSTs extracted from maize participated in herbicide metabolism by facilitating the conjugation of synthetic chemicals with glutathione [GSH; γ-Glu–Cys–Gly] [2].Plant GST proteins are globular dimers, typically consisting of two separate N- and C-terminal domains [3]. The N-terminal domain has a structural configuration of β-strands and α-helices organized in a thioredoxin-like fold, with a GSH binding site (G-site) characterized by excellent substrate selectivity. The C-terminal domain is composed solely of α-helices, including an electrophilic substrate binding site (H-site) characterized by hydrophobic properties and extensive substrate specificities. This region accommodates many substrates and ligands [1,4]. Previous studies have classified plant GSTs categorizes into 14 classes, including phi, zeta, tau, theta, lambda, Ure2p, EF1Bγ, hemerythrin, iota, glutathionyl-hydroquinone reductases (GHRs) [5], dehydroascorbate reductase (DHAR), tetrachlorohydroquinone dehalogenase (TCHQD), metaxin, and microsomal prostaglandin E synthase type 2 (mPGES-2) [6,7,8]. Among these, tau, phi, DHAR, and lambda are plant-specific GST classes. The tau and phi classes are the most prevalent, likely due to their rapid development in plants. The Arabidopsis thaliana GST gene family consists of 57 members, categorized into nine classes, including 13 phiand 28 tau members [9]. The functions of GST classes apart from tau and phi remain poorly understood. DHAR GSTs diminished the production of ascorbic acid within the ascorbate–glutathione cycle to produce antioxidants [10]. GSTs of the zeta class exhibit isomerase activity and are involved in tyrosine catabolism [11]. Theta class GSTs function as glutathione peroxidases, contributing to the metabolism of oxidative stress [12]. Lambda GSTs use flavonols as high-affinity ligands. Dixon and Edwards [4] elucidated a unique glutathione-dependent function of these enzymes in the recycling of oxidized quercetin. DHAR and lambda GSTs exhibited action towards dehydroascorbic acid (DHA) [13].

Previous studies have also demonstrated the functional diversity among GST genes [14,15]. Despite limited understanding, iota class GSTs appear to function as oxidoreductases rather than conjugating GSTs, while hemerythrins facilitate the detoxification of heavy metals by promoting the conjugation of GSH with metal ions. Recent research has demonstrated the scavenging ability of members of tau class against most hazardous reactive carbonyl species (RCS) (i.e. acrolein) and proposed their participation in defense mechanism against RCS [16].GSTs primarily facilitate the conjugation of reduced tripeptide glutathione (GSH) to diverse reactive electrophiles. GST proteins are crucial for the appropriate growth and physiological functioning of plants. They are crucial in stress response [17,18] and participate in the breakdown of insecticides and herbicides [18]. GSTs exhibit responsiveness to plant hormones, including auxins, ethylene, salicylic acid, abscisic acid, and jasmonic acid [10,19]. Based on current advances in plant proteomics, genomes, and transcriptomics studies, Chronopoulou et al. (2017) examined the functions of GSTs. Apart from their catalytic activity, GSTs bind to flavonoids and anthocyanins, therefore acting as non-catalytic proteins by moving them from the cytoplasm into the central vacuole [5].Involved in plant development, GSTs have binding affinity to auxin [20] and cytokinin [21], phytohormones. Additionally binding to porphyrinogens [22] and oxylipins [9], GSTs guard plant cells against oxidative stress [23,24]. Furthermore crucial in responses to biotic and abiotic stressors like pathogen assault, heavy metals [25], dehydration [26] and salt [27], as well as salicylic acid signalling [13] are plant GSTs. Still, the specific processes are unknown.

Furthermore, GSTs may play a role in the transport and metabolism of secondary metabolites, including anthocyanins, flavonoids, and porphyrins [4]. Recent genome-wide analyses have identified 53 GST genes in Arabidopsis [24], 79 in rice [14,28], 25 in soybean, 42 in maize [29], 90 in tomato [30], 85 in pepper [31], 32 in pumpkin [32], 23 in sweet orange [33], 90 in potato [34], and 20 in Dracaena cambodiana [35]. Despite the identification of Glutathione S-Transferases (GSTs) as contributors to plant stress tolerance and detoxification, a comprehensive analysis of the GST gene family in Crocus sativus (saffron) remains unavailable. The absence of genetic and functional insights limits our understanding of the roles saffron GSTs play in stress adaption and overall plant resilience. Understanding GST operations may facilitate improvements in saffron production and quality among environmental challenges. This study conducted a genome-wide annotation and phylogenetic analysis of the GST gene family in Crocus sativus. Using computational tools such as BLAST [36], TBtools [37], and Python scripts, we identified and characterized GST genes, examined their motif structures, and assessed their physicochemical characteristics. Comparative analyses with Arabidopsis thaliana and other plant species revealed lineage-specific gene divergence and evolutionary trends. Promoter region analysis and phylogenetic tree construction provided deeper insights into the regulatory and functional diversity of GST genes. This study enhances our understanding of the evolution and functional divergence of GST gene families in Crocus sativus and other monocots, offering a foundation for future research on plant stress tolerance and metabolic adaptation.

2. Materials and Methods

2.1. Plant Materials

The study "Crocus genome reveals the evolutionary origin of crocin biosynthesis" [38] provided the genomic dataset of Crocus sativus used in this work. This dataset, which includes annotated genes, protein sequences, and genomic scaffolds, offers a comprehensive toolkit for identifying and studying the function roles of GST genes. Using GST protein sequences of Arabidopsis thaliana, we retrieved reference sequences, including NP_173160.1, AEE29554.1, and AAD34992, from NCBI databases for a comparative analysis. These sequences were curial for identifying homologous GSTs in Crocus sativus and evolutionary research. Building on a solid foundation created by merging the Crocus sativus genome with publicly available data from Arabidopsis thaliana, we began investigating the lineage-specific evolutionary and functional variations within GST gene families.

2.2. Identification and Annotation of GST Genes

Using known GST protein sequences from Arabidopsis thaliana (e.g., NP_173160.1, AEE29554.1, and AAD34992.1) as query sequences, BLASTp searches were performed against the annotated protein sequences of the Crocus sativus genome dataset to identify GST genes. High-confidence candidates were selected based on a stringent e-value cutoff (<1e-5) and the presence of GST-specific conserved domains, as confirmed by HMM-based searches using the Pfam [39] and SMART databases [40]. Redundant sequences were removed, and annotations of the selected GST genes were validated using TBtools [37]. The identified protein sequences were compiled into a FASTA [41] file for downstream analyses, including motif characterization, phylogenetic analysis, and structural studies. These procedures ensured the construction of a high-quality GST gene for Crocus sativus collection. Additionally, GST genes from Arabidopsis thaliana and other allied plant species were analyzed using the same pipeline to establish an evolutionary framework.

2.3. Phylogenetic Analysis

The evolutionary relationships of GST proteins in Crocus sativus were investigated using phylogenetic analysis. ClustalW [42] algorithm in MEGA [43] matched protein sequences of GST genes found in Crocus sativus with orthologous GST protein sequences from Arabidopsis thaliana and other related species. Multiple sequence alignment guaranteed exact comparison between conserved and divergent areas among the sequences.A phylogenetic tree was constructed using the Neighbor-Joining (NJ) [44]method in the MEGA program. The reliability of tree branches was assessed using 1,000 bootstrap replicates, with default settings including pairwise deletion for gaps. GST proteins were classified into several subfamilies (e.g., Tau, Phi, Zeta), and evolutionary events such as gene duplication and divergence were inferred from the resultant tree. This study provides a framework for understanding the functional specialization and lineage-specific evolution of the GST gene family in Crocus sativus.

2.4. Motif Analysis

The MEME Suite [45] was used to find conserved motifs among GST proteins of Crocus sativus. Using default settings for motif width and significance thresholds, the protein sequences found in the previous phase were analyzed under criteria allowing the detection of up to 15 motifs. By cross-referencing these motifs with domain annotations, their functional importance was evaluated.

The functional domains of GST proteins were confirmed using the Pfam and SMART databases, which identified characteristic GST domains, including C-terminal alpha-helical domain and the N-terminal thioredoxin-like domain. The processed data were systematically organized to facilitate visualization and further interpretation in the results stage. These analyses provided insights into the conserved structural elements and functional roles of GST proteins in Crocus sativus.

2.5. Domain Validation

Domain validation using the Pfam database was performed to confirm the identification of GST genes in Crocus sativus. Conserved domains characteristic of the GST family were identified, including the alpha-helical C-terminal domain (GST_C) and the thioredoxin-like N-terminal domain (GST_N). Reflecting the diversity of GST proteins in Crocus sativus, variations of these domains, such as GST_N_2, GST_N_3, and GST_N_4, as well as GST_C_2, GST_C_3, GST_C_7, and GST_C_6 were also found. Additionally, certain sequences contained other domains associated with functional specialization, such as DnaB_bind and CLIC-like_N. These domain variations highlight the structural and functional diversity of the GST gene family in Crocus sativus. Only genes containing conserved GST-specific domains were included in the final dataset for further phylogenetic and functional analysis.

2.6. Ka/Ks Analysis of GST Gene Sequences

Synonymous (Ks) and non-synonymous (Ka) substitution rates were calculated to evaluate the selective pressures acting on the GST gene family in Crocus sativus. To ensure high-quality alignment for subsequent analysis, coding sequences (CDS) of the identified GST genes were aligned using the ClustalW method in the MEGA program. Subsequently, a neighbor-joining phylogenetic tree was constructed from the aligned sequences to establish an evolutionary framework and guide pairwise comparisons.

The Codon-Based Test of Selection in MEGA computed pairwise Ka/Ks values. The Ka/Ks ratio of each gene pair provided insight into the type of selection acting on the sequences. Ratios greater than 1 indicated positive selection, suggesting adaptation to specific functional or environmental pressures. Ratios less than 1 suggested purifying selection, which works to preserve gene functionality by removing harmful mutations. A ratio equal to 1 indicated neutral evolution. Finally, the pairwise distance data, including the Ka/Ks ratios, were saved in an Excel file for further analysis. This approach facilitated the assessment of the evolutionary processes influencing the GST gene family in Crocus sativus as a whole.

2.7. Location Analysis of Chromosomes

The genomic structure of GST genes in Crocus sativus was analyzed through chromosomal distribution mapping. Genomic annotation data, including the locations of GST genes, was retrieved from the Crocus sativus genome. Using TBtools, the GST genes were mapped to the chromosomes. The mapping process involved charting the physical positions of GST genes on each chromosome and converting the genomic data into an appropriate format. A phenogram was created by grouping GST genes based on their chromosomal proximity or clustering patterns, thus revealing their genetic architecture. Each group was color-coded to enhance the visualization of their distribution across chromosomes. This analysis identified gene clusters that could be evolutionary hotspots or regions affected by gene duplication events. The generated graphics clearly illustrated the spatial arrangement and positions of GST genes within the Crocus sativus genome.

2.8. Predicting Subcellular Localization

WolfPsort [46], a computational technique designed to predict protein localization based on amino acid sequences, was used to project the subcellular localization of GST proteins in Crocus sativus. The amino acid sequences of all 60 identified GST genes were submitted to the WolfPsort server. The tool analyzed the sequences and assigned a score for each protein across various cellular compartments, including the cytoplasm, chloroplast, nucleus, vacuole, mitochondria, cytoskeleton, extracellular space, and others. These scores represent the likelihood of a protein localizing to a particular compartment. The compartment with the highest score was considered the primary localization site for each protein. If significant scores were observed, secondary and tertiary localization predictions were also considered to provide a comprehensive view of the subcellular location of each GST protein. The resulting data were compiled for further display and analysis.

2.9. Physicochemical Property Analysis

The structural and functional characteristics of the identified GST proteins in Crocus sativus were analyzed through physicochemical property assessments. The ExPASy [47] ProtParam tool was used to compute key parameters, including molecular weight (kDa), isoelectric point (pI), and hydrophobicity (GRAVY score). FASTA-format protein sequences were input into ProtParam, and the results for each sequence were recorded. The GST proteins were classified based on their physicochemical profiles, including their acidic or basic nature (based on pI) and their hydrophilic or hydrophobic tendencies (based on the GRAVY score). This analysis is crucial for predicting the biochemical functions of GST proteins and selecting candidates for further functional studies. By providing insights into the stability, solubility, and functional significance of these proteins, this analysis lays the groundwork for experimental validation and comparative studies.

2.10. Aligning

Multiple sequence alignments using TBtools were performed to examine the sequence conservation among GST genes. FASTA format protein sequences were aligned to identify conserved regions and motifs. Color schemes in the alignment facilitated visualization by highlighting both similarities and differences in amino acid residues. Accurate representation was ensured using the default settings, and the data were exported for further annotation and analysis.

2.11. Hormone Response Element Promoter Analysis

Using genome annotation data, a 1500 bp upstream sequence from the transcription start site (TSS) was extracted to identify hormone-responsive cis-regulatory elements in the promoter regions of the discovered Crocus sativus GST genes using PlantCare website [48]. The TBtools software package was used for promoter region analysis. Specifically, hormone-responsive motifs such as MEJA (Methyl Jasmonate), ABRE (Abscisic Acid Responsive Element), and GARE (Gibberellin Responsive Element) were selected for study. Pattern matching against a known database of hormone-regulating sequences formed the basis for cis-element identification. To visualize their presence and distribution throughout the promoter regions of each gene, graphical representations of the hormone-responsive element density were created.

2.12. Protein Modeling and Functional Analysis Using Swiss-Model and PDBSum Tools for GST

The target protein, glutathione S-transferase (GST), has an amino acid sequence obtained from a genomic database. The sequence was analyzed using the Swiss-Model [49] Interactive Workspace for structural modeling. After submitting the sequence for template search, the Swiss-Model server identified homologous templates with high structural similarity. Further analysis focused on the template with the highest sequence identity and query coverage. The Tau class GST from Mangifera indica, with a 56.89% sequence identity and 0.96 coverage (PDB ID: 5g5e), was selected as the template.X-ray diffraction at a 2.35 Å resolution was used to solve the structure, providing detailed information on ligand binding and the active site. PDBSum [50] was used to identify the N-terminal and C-terminal domains, thus confirming the structural features of the protein. Using the Swiss-Model and PDBSum tools, the selected model was further investigated for ligand-binding interactions and functional domains. Key ligands, including S-hexylglutathione and di(hydroxyethyl)ether, were examined for interactions with active site residues within a 4 Å radius.

2.13. Glutathione S-Transferase (GST) Gene Expression Variations in Source and Sink Tissues During the Development of Saffron (Crocus sativus L.) Corms

Raw sequencing data for this study were obtained from SRA run accessions SRR24751341 toSRR24751364, while RNA sequencing (RNA-Seq) data retrieved from the NCBI Sequence Read Archive (SRA) under project accession PRJNA976833.The datasets were generated from saffron (Crocus sativus L.) samples collected at various developmental stages, including both source (leaf) and sink (corm) tissues. Raw sequencing data were analyzed and their quality assessed using FastQC (v0.11.9). Low-quality reads and adapter sequences were trimmed using Trimmomatic (v0.39), ensuring high-quality sequences for further analysis. HISAT2 (v2.2.1) was used to align the processed reads to a custom saffron transcriptome, and SAMtools (v1.11) [51] was applied to convert the resulting SAM files into BAM format. FeatureCounts (v1.6.5) [52] was used to quantify gene expression levels, and DESeq2 (v1.30.1) in R was employed for differential gene expression (DEG) analysis. Genes with an adjusted p-value < 0.05 and a log2 fold change > 2 were considered significant. The expression patterns of glutathione S-transferase (GST) genes across various tissues and developmental stages were visualized using a heatmap.

3. Results

3.1. Discovery of GST Genes in Crocus sativus

A comprehensive genome-wide study identified 60 GST genes in Crocus sativus. These genes are distributed across the genome and annotated based on the presence of conserved GST-specific domains, such as the C-terminal alpha-helical domain and the thioredoxin-like N-terminal domain. Comprehensive genomic data for the identified Crocus sativus GST genes are provided in Supplementary Table S1, which includes key information such as gene ID, gene name, chromosome, strand orientation (start and end), and amino acid count. This dataset offers a detailed overview of the genes, highlighting their genomic locations, structural features, and potential functions, serving as a foundational resource for further research on the GST gene family in Crocus sativus.

3.2. Phylogenetic Analysis of GST Genes

A comprehensive phylogenetic analysis was conducted to explore the functional diversity and evolutionary relationships of GST genes in Crocus sativus. Protein sequences of the identified Crocus sativus GST genes were selected based on their functional relevance in previous studies and aligned with well-characterized GST genes from other plant species. These included CkmGST3 from Cyclamen persicum (associated with anthocyanin storage), AtGSTU19 and AtGSTF10 from Arabidopsis thaliana (linked to abiotic stress tolerance and oxidative stress response), and LcGST4 from Litchi chinensis (involved in anthocyanin transport to vacuoles for pigmentation). These genes were retrieved from NCBI using their respective accession numbers: LcGST4 (ALY05897.1), AtGSTU19 (Q9ZRW8.1). Their roles in various plant physiological processes are well-documented in previous research [53,54].

Motif and domain analysis confirmed the presence of GST-specific conserved domains, including GST_N, GST_C, and their variants (e.g., GST_N_3, GST_C_7), highlighting their functional importance. Sequence alignment was performed using MEGA, and a neighbor-joining phylogenetic tree was constructed. TBtools was used for tree visualization. Specific clusters were identified that linked Crocus sativus GSTs with reference genes. For example, CSGST-14 grouped closely with LcGST4, suggesting a conserved role in anthocyanin transport, while CSGST-21 and CSGST-35 clustered with AtGSTU19, indicating potential involvement in detoxification and stress response pathways. This study provided insights into the evolutionary pathways and functional diversity of GST genes in Crocus sativus. Although indicating unique adaptations in Crocus sativus GST genes, the inclusion of reference genes from other species demonstrated shared functional characteristics such stress tolerance and secondary metabolite transport. Supplementary Table S2 includes a list of glutathione S-transferase (GST) genes discovered in Crocus sativus, along with their gene identifiers, protein family classifications (Pfam), gene locations with statistical significance (E-value), and brief descriptions of their functions. These genes belong primarily to the glutathione S-transferase family, which is involved in cellular detoxification processes. They are categorized based on the enzymes they encode, such as TAU, PHI, THETA, and others. This dataset provides critical information regarding the structure and classification of GST genes in Crocus sativus.

3.3. Motif Interpretive Analysis of GST Genes

A detailed investigation of the GST genes in Crocus sativus revealed conserved sequence patterns critical for their structural and functional roles. Using the MEME Suite, fifteen distinct motifs were identified, and their positions along with corresponding Pfam values were recorded. These motifs, crucial for maintaining GST-specific enzymatic activity and binding properties, highlight conserved regions common across GST family members. Every known GST gene contains motifs, some of which are ubiquitous, while others are restricted to specific gene sets. Supplementary Table S3 presents the sequence and width of each motif.

3.4. Domain Analysis Conservation in GST Genes

Motif and domain analysis of the identified Crocus sativus GST genes revealed several conserved domains unique to the GST family. Consistent with canonical GST activity and structural stability, the N-terminal GST domains (e.g., GST_N, GST_N_2, GST_N_3, GST_N_4) and C-terminal GST domains (e.g., GST_C, GST_C_2, GST_C_3, GST_C_6) were predominant throughout the sequences. In addition to these core GST domains, some sequences exhibited less common or more specialized domains, such as CLIC-like_N, DnaB_bind, and EF1G, potentially indicating functional divergence or specialization within the GST family.

As shown in Supplemental Table S2, the genes were further categorized based on their domain composition and phylogenetic relationships. The THETA group contains genes with GST_C and associated C-terminal domains, while the TAU group consists of genes with GST_N and GST_N-3 domains. Genes with less common domains, such as GST_N_2 and GST_N_4, are classified under the PHI category, and those with unique or unknown domains fall into the OTHER group. These classifications provide insight into the functional and evolutionary aspects of the GST gene family.

Figure 1 focuses on motif and domain structure. It depicts the genomic architecture of various gene sequences, including CSGST and GST family members. The evolutionary connections among different isoforms (e.g., CSGST-30-Theta, CSGST-31-Theta) are represented in a tree-like structure, showing their association with specific motifs or protein domains. The motifs are color-coded based on protein families, including MS_channel, glutaredoxin, atrophin-1, and the GST superfamily. These motifs are mapped along the gene sequences, with each bar representing a different domain within the sequences. The image on the right shows additional structural annotations, including coding sequences (CDS) with conserved motifs and potential functional regions across the isoforms. The variability in motif distribution observed in this study may reflect functional specialization or evolutionary divergence within the gene family. Among the significant protein domains relevant to functional studies are the GST and Atrophin-1 superfamilies, which could guide further research on gene functions in adaptive mechanisms or biochemical processes.

In Figure 2, Panel a presents an evolutionary tree illustrating the relationships among GST protein sequences from various species and isoforms (e.g., CGST and ATS). The tree highlights the phylogenetic groupings, with color-coding of different protein clades. This visualization underscores the structural and functional differences across protein variants, emphasizing the evolutionary diversification of GST proteins. Panel b shows sequence logos for the C-terminal (top) and N-terminal (bottom) regions of GST proteins. These logos encode the amino acid distribution at each position within the corresponding domains, with the height of each character representing the relative frequency of the matching amino acid. This provides a quantitative view of sequence conservation and variation. The visualization highlights the functional relevance and evolutionary conservation of the terminal regions, which are crucial for GST catalytic activity and structural integrity.

Taken together, these visualizations offer a comprehensive comparative study of GST protein evolution and sequence architecture, providing valuable insights into the functional implications of sequence variation and domain-specific conservation across species.

3.5. Ka/Ks Analysis of GST Gene Sequences

The Ka/Ks analysis of GST genes in Crocus sativus provided valuable insights into the evolutionary pressures shaping this gene family. Most gene pairs exhibited Ka/Ks ratios smaller than 1, indicating strong purifying selection that maintains the functional integrity of these genes. This suggests that most GST genes have evolved under constraints to preserve their essential roles in stress response and detoxification pathways. Gene pairs with Ka/Ks ratios greater than 1 indicated positive selection, suggesting adaptive changes to specific environmental or metabolic conditions. These gene pairs likely contribute to the functional diversity of the GST family by engaging in specialized functions or evolutionary divergence. Supplementary Table S4 presents the pairwise Ka/Ks values, highlighting the predominance of purifying selection and the occasional occurrence of positive selection within the Crocus sativus GST gene family. The table provides an evolutionary analysis of gene pairs, focusing on synonymous (KS) and nonsynonymous (KA) changes. The Ka/Ks ratio is used to assess the evolutionary pressure on these genes: values less than one indicate purifying (negative) selection, values around one suggest neutral evolution, and values greater than one indicate positive (diversifying) selection. Each gene pair is categorized according to its evolutionary selection type, such as "Negative or Purifying" or "Positive or Diversifying," to explain the origins of these genes within the species.

3.6. GST Genes Chromosomal Distribution

Figure 3 illustrates the distribution of GST genes across the eight chromosomes of Crocus sativus. While some chromosomes, particularly Chr1 and Chr2, contained a higher number of GST genes, others had fewer. On several chromosomes, genes were concentrated in specific regions, with notable hotspots observed on Chr8. These gene clusters likely reflect the results of gene duplication events, which may have facilitated the expansion and diversification of the GST gene family in Crocus sativus. The distribution pattern, with separate groups color-coded in the phenogram, also highlights the functional organization of GST genes. Group 1 (red), Group 2 (blue and purple), and Group 3 (green) suggest a potential functional or structural classification of the GST genes. These findings offer valuable insights into the evolutionary processes and genomic architecture of the GST gene family in Crocus sativus.

3.7. GST Proteins Subcellular Localization

Thesubcellular localization of GST proteins in Crocus sativus was analyzed using the WolfPsort tool, revealing diverse distributions across various cellular compartments, reflecting their multiple roles within the plant. Most GST proteins, including CSGST-1, CSGST-2, and CSGST-16, were predicted to localize primarily in the cytoplasm, supporting their involvement in detoxification and cellular redox maintenance.

Several GST proteins, such as CSGST-15 and CSGST-20, were predicted to localize in the vacuole, suggesting roles in the transport or sequestration of secondary metabolites. Proteins like CSGST-27 and CSGST-29 showed localization to the nucleus and mitochondria, implying potential regulatory functions or involvement in stress responses. Additionally, proteins such as CSGST-13 and CSGST-17 were predicted to be localized in the chloroplast, possibly linked to photosynthesis.

Supplementary Table S5 provides a detailed overview of the subcellular distribution of the CSGST genes, including their localization in the cytoplasm, nucleus, chloroplast, and mitochondria. The localization patterns of these proteins highlight their diverse functional roles in cellular metabolism and stress response pathways, with some proteins showing distribution across multiple compartments.

3.8. Physicochemical Characteristics of GST Proteins

The physicochemical properties of the GST proteins in Crocus sativus exhibited a range of structural and functional characteristics. The molecular weights of these GST proteins ranged from 20.49 kDa to 68.49 kDa, while the isoelectric points (pI) varied between 5.15 and 7.56, indicating a predominance of neutral to moderately acidic proteins. The Grand Average of Hydropathy (GRAVY) scores ranged from -0.281 to 0.143, suggesting that most GST proteins are hydrophilic, consistent with their enzymatic activity in aqueous environments. Based on subcellular localization predictions, most GST proteins were identified in the cytoplasm, with several showing dual localization to both the cytoplasm and chloroplast. These results highlight the functional diversity of GST proteins in Crocus sativus and suggest their potential roles in maintaining redox balance, responding to stress, and facilitating cellular detoxification processes. Supplementary Table S6 provides a comprehensive list of the molecular characteristics and subcellular localization of CSGST proteins, including Gene ID, protein molecular weight (kDa), isoelectric point (pI), GRAVY score, and predicted subcellular localization. While most of the CSGST proteins are expected to be cytoplasmic, others are predicted to localize in the chloroplast, mitochondria, nucleus, or plasma membrane, reflecting their specialized roles in cellular processes. Most of these proteins exhibit hydrophilic properties based on the GRAVY scores.

3.9. Alignment

The multiple sequence alignment (MSA) of GST proteins revealed highly conserved regions that correspond to their functional domains. Notably, conserved regions matched the N-terminal and C-terminal GST domains, highlighting their evolutionary importance in maintaining catalytic activity and structural integrity. Variations in non-conserved regions suggest functional diversity among the GST members. The alignment also highlighted unusual residue patterns in several sequences, which may indicate adaptations to specific biological activities or stress responses. The alignment is color-coded to emphasize conserved motifs and variants. Figure 4 presents the multiple sequence alignment (MSA) of protein sequences derived from CSGST genes. Various colors highlight conserved amino acid residues, with the letter height in the sequence logo reflecting the degree of conservation. The sequence logo underscores the critical conserved sections that are essential for protein activity, illustrating the frequency of specific amino acids at each position. This work provides a detailed analysis of structural elements, sequence conservation, and functional domains within the aligned sequences.

3.10. Hormone-Responsive Element Distribution in Promoter Regions

Analysis of the promoter regions of Crocus sativus GST genes revealed the presence of hormone-responsive cis-elements, including MEJA, ABRE, and GARE motifs, with varying distribution patterns. MEJA elements, which are linked to methyl jasmonate regulation, were the most abundant and prominently located in the upstream sequences of genes such as Csativus01200 (CSGST-14), Csativus27729 (CSGST-21), and Csativus51083 (CSGST-7), as illustrated in the distribution map (Figure 5). ABRE motifs, associated with abscisic acid response, were predominantly identified in genes like Csativus40095 (CSGST-26) and Csativus55558 (CSGST-52). Although less common, GARE motifs were also found in genes such as Csativus57966 (CSGST-16) and Csativus44106 (CSGST-27), indicating potential involvement in gibberellin-mediated regulation. The positioning and abundance of these motifs suggest their crucial role in regulating GST gene expression, contributing to plant stress responses and growth through complex hormonal signaling networks.

3.11. Protein Modeling and Functional Analysis of Glutathione S-Transferase Using Swiss-Model and PDBSum Tools

The structural modeling of the GST protein from Crocus sativus (sequence: MSNNGAGEEVKLLGLWASPYVIRARIALNIKGVDHEFLEETFGRKSDLLIESNPVYKKMPVLIHQGKTVCESMIIVQYIDDVWNSAKPSIVPSDPYELAMARFWAVYIDDNWTSSVIGMMRATTEEARVAAIEQAAKALQTLENAFEKHYKGKLFFGGEKIGYLDIALGSSVGWIGVLEKISGVTLLDEARTPQLMGWVARFWSQDVVRSVIPETDRFVEFSKEVFATRAAQAAK, Sequence length: 235 residues) was conducted with the template PDB ID 5g5e, corresponding to a Tau-class GST from Mangifera indica. This template, with 56.89% sequence identity and 96% query coverage, proved to be a reliable foundation for accurate structural representation. Structural analysis revealed two critical domains: an N-terminal domain facilitating substrate binding and a C-terminal domain regulating catalytic activity. These domains were visualized using the Swiss-Model platform and further validated with PDBSum annotations. Functional site analysis showed interactions between the GST protein and the ligand S-hexylglutathione, including 11 residues within a 4 Å radius that formed hydrogen bonds and hydrophobic interactions. Non-functional ligand di(hydroxyethyl) ether suggested additional binding sites in the active pocket. The modeled GST protein structure is homo-dimeric, with a QMEANDisCo global score of 0.80 ± 0.05, indicating excellent structural reliability. These findings provide valuable insights into the structural and functional characteristics of the GST protein. Figure 6 illustrates the structural and functional aspects of the Tau-class glutathione S-transferase (GST) from Mangifera indica, including: (a) a 3D structural representation, (b) the enzyme's role in detoxification via glutathione conjugation, (c) key conserved domains essential for GST function, (d) the secondary structure composed of β-sheets and α-helices, (e) the importance of two distinct structural domains in enzymatic activity, and (f) its involvement in redox homeostasis and detoxification pathways, as shown in the KEGG pathway for glutathione metabolism.

3.12. Glutathione S-Transferase (GST) Gene Expression Variations in Source and Sink Tissues During the Development of Saffron (Crocus sativus L.) Corms

The analysis of RNA-Seq data revealed significant differential expression of GST genes between source (leaf) and sink (corm) tissues in saffron. A subset of GST genes was notably upregulated in corm tissues, particularly during the bulking and maturation phases, indicating their involvement in key metabolic and stress response pathways critical for corm development. In contrast, GST gene expression in leaf tissues exhibited varied patterns, with some genes showing increased expression during the early stages of corm development. The heatmap visualization in Figure 7 illustrates the differential expression levels of GST genes, with red hues representing upregulated genes, and expression values ranging from -3.00 to +5.00, reflecting the extent of gene regulation across different samples. These findings underscore the essential role of GST genes in the molecular regulation of saffron corm growth, particularly in relation to environmental stress and metabolic demands.

4. Discussion

The findings of this study provide a comprehensive understanding of the GST gene family in Crocus sativus, offering significant insights into their role in plant stress responses and secondary metabolism. The identification of 60 GST genes with conserved GST_N and GST_C domains supports the hypothesis that these genes are crucial for secondary metabolites transport and detoxification, functions well-documented in other plant species such as Arabidopsis thaliana and Vitis vinifera [55]. The observed gene duplications align with previous studies, suggesting that these duplications contribute to the functional diversity of GSTs in response to environmental stresses [19]. Additionally, the presence of hormone-responsive motifs in the promoters of GST genes indicates their potential involvement in mediating hormone-regulated stress responses. This role is further supported by studies of other plants susceptible to abiotic stresses, such as oxidative damage and drought. The subcellular localization analysis confirms the involvement of many GSTs in cellular detoxification systems, with several proteins detected in the cytoplasm, chloroplasts, and vacuoles. Based on their localization patterns, which suggest various plant functions, it is hypothesized that GSTs may play roles beyond stress responses, including in the transport of secondary metabolites and other cellular processes.

The localization pattern of GSTs, along with their diverse functions in plants, supports the hypothesis that these enzymes participate in cellular processes beyond stress regulation, including the synthesis and transport of secondary metabolites. Phylogenetic analysis of the GST family revealed both conserved and species-specific evolutionary patterns, suggesting that while certain functions are conserved across plant species, the Crocussativus GSTs may have evolved specific roles contributing to the plant's unique biochemical pathways, such as the biosynthesis of crocin. From an agricultural perspective, these findings provide valuable insights into the molecular mechanisms that could be harnessed to enhance saffron's stress tolerance and productivity, essential for sustaining its high-value cultivation amid changing climatic conditions. The identification of conserved functional domains involved in stress responses and detoxification pathways opens avenues for future research aimed at genetically improving saffron’s resilience. Functional validation of these GSTs through gene-editing technologies like CRISPR, along with investigations into their interactions with other stress-related proteins, should be a priority for future studies. Furthermore, exploring the role of GSTs in the biosynthesis of crocin and other bioactive compounds could help optimize saffron cultivation and maximize the production of these valuable metabolites.

5. Conclusions

This comprehensive study focused on the identification, annotation, and functional analysis of glutathione S-transferase (GST) genes in Crocus sativus. Sixty GST genes were identified through genome-wide searches and categorized into distinct subfamilies based on conserved domains, motif characteristics, and evolutionary relationships. Chromosomal distribution analysis revealed clustering of GST genes, indicating gene duplication events. Ka/Ks analysis suggested that purifying selection is the primary evolutionary driver, with sporadic positive selection pointing to functional diversity. Promoter analysis uncovered hormone-responsive cis-elements, including MEJA, ABRE, and GARE, highlighting the role of GST genes in hormone-mediated stress responses. Predictions regarding subcellular localization showed that GST proteins are predominantly localized in the cytoplasm, with occasional presence in the nucleus, vacuole, or chloroplast, emphasizing their versatile functions. The structural characteristics of a GST protein were analyzed using the Swiss-Model platform and PDBSum tools, revealing the N-terminal and C-terminal domains responsible for substrate binding and catalytic activity. The reliability of the model was confirmed through ligand interaction analysis, which identified key residues involved in functional connections. These findings enhance our understanding of the structure, evolution, and function of GST genes in Crocus sativus, and open new avenues for research on GST-mediated stress tolerance, plant growth, and the enhancement of crop resilience and secondary metabolite synthesis.