Glutathione S-transferase in mediating adaptive responses of oats (Avena sativa) to osmotic and cadmium stress: genome-wide analysis

Abstract

Background

Glutathione S-transferases (GSTs) are essential multifunctional enzymes. In the face of abiotic stresses such as drought and heavy metal exposure, plants utilize GSTs for detoxification and antioxidant defense, as these enzymes facilitate the conjugation of glutathione (GSH) with toxic compounds. Specific details of this process, however, remain unknown.

Results

This study identified 118 Avena sativa GST (AsGST) genes within the A. sativa genome and classified them into five subfamilies: Tau, Phi, Zeta, Lambda, and EF1Bγ. Phylogenetic analysis revealed that AsGSTs exhibit significant similarity to corresponding GST categories in Arabidopsis thaliana and Oryza sativa, indicating a possible common ancestor. Gene structure and conserved motif analysis demonstrated that AsGST genes within the same subfamily shares similarities in the number and positioning of exons and introns, as well as in motif composition, suggesting that these genes may perform analogous biological functions in A. sativa. The promoter regions of the identified genes are enriched with various cis-acting elements that play roles in plant growth and development, stress response, and hormone signaling. Transcriptomic analysis and real-time quantitative PCR (RT-qPCR) validation indicated that the expression of four AsGST genes (AsGSTU12, AsGSTU13, AsGSTU14, and AsGSTU15) was significantly up-regulated in the roots of A. sativa under both PEG-induced drought stress and CdCl₂-induced cadmium stress. These genes likely regulate reactive oxygen species (ROS) levels by catalyzing their scavenging through glutathione (GSH) substrates, and may also participate in ABA signaling and the maintenance of osmotic homeostasis. Under cadmium stress, these genes may mitigate cadmium toxicity by enhancing the chelation and sequestration of cadmium via GSH or through its compartmentalization, as evident from the subcellular localization studies.

Conclusion

This study systematically described the GST gene family in A. sativa, characterized its expression patterns and potential functions in response to drought and cadmium stress, and confirmed the essential role of the AsGST gene family in mediating stress responses. The findings enhance our understanding of the mechanisms underlying stress tolerance and offer valuable genetic resources for breeding stress-tolerant A. sativa. The work also provides a theoretical framework and identifies gene targets for the development of stress-resistant A. sativa varieties.

Background

Glutathione S-transferases (GSTs) represent a diverse group of multifunctional enzymes prevalent in plants, where they are essential for growth, development, and responses to environmental stresses [1]. These enzymes facilitate various biological processes, including detoxification metabolism, antioxidant defense, and signaling, by catalyzing the conjugation of glutathione (GSH) to electrophilic compounds [2,3,4]. A typical GST comprises two distinct binding sites, which are located at the N- and C- terminus, to bind GST (GST-N) and neighboring substrate (GST-C), respectively [5, 6]. Based on similarities in the sequence and structural as well as their functional roles, plant GST genes are categorized into 14 subfamilies, including Phi (F-type), Tau (U-type), Lambda (L-type), Theta (T-type), Zeta (Z-type), dehydroascorbate reductase (DHAR), γ-subunit of the eukaryotic translation elongation factor 1B (EF1Bγ), tetrachlorohydroquinone dehalogenase (TCHQD), metaxin, Ure2p, Hemerythrin (H), Iota (I), microsomal prostaglandin E-synthase type 2 (mPGES- 2), and glutathionyl-hydroquinone reductase (GHR) [7]. Notably, the Tau, Phi, Lambda, and TCHQD subfamilies are exclusive to plants and play pivotal roles in xenobiotic metabolism by binding a broad spectrum of harmful exogenous compounds, such as insecticides and herbicides, thus safeguarding plants from toxic effects [8].

The multifunctionality of GSTs is evidenced by their involvement in regulation plant ROS signalling, redox balance, and hormone biosynthesis [9]. GSTs are intricately linked to biotic and abiotic stresses [10] and are vital for plant growth, development, and enhanced stress tolerance, playing a crucial role in the antioxidant defense mechanism [11]. The expression and activity of GSTs significantly bolster plant adaptation to various stress conditions, including low temperatures, high salinity, drought, and heavy metal exposure [12]. For instance, under salt stress, GSTs in Gossypium hirsutum exhibit tissue-specific expression patterns [13], while the SlGSTU43 gene in Solanum lycopersicum aids in scavenging ROS and promoting lignin biosynthesis, both of which enhance salt stress tolerance [14]. In poplar, PtGSTF1 improves ion homeostasis and ROS scavenging under salinity stress [15] while in Juglans regia the JrGSTTau1 gene enhances cold tolerance [16]. Similarly, the IbGST4 and IbGST2 genes mitigate low-temperature-induced ROS accumulation and associated damage [17].

Also significant is the role of GST in responding to drought and heavy metal stress. For example, CsGSTU8 in Camellia sinensis positively regulates drought stress, modulated by the transcription factor (TF) CsWRKY48 [18]. In Pisum sativum GST activity was elevated in plants subjected to cadmium stress [19]. Ectopic expression of OsGSTU4 and OsGSTU30 from rice has been shown to enhance the tolerance of transgenic Arabidopsis thaliana to salinity, oxidative stress, drought, and heavy metal stress [20, 21]. These findings highlight the growing importance of GST in addressing drought and heavy metal stress, suggesting potential strategies for the development of resilient crops capable of thriving in challenging environmental conditions. However, given the plethora of GST isoforms, it remains to be answered what are their specific roles, and which of them may be a suitable target for improving abiotic stress tolerance in crops.

Avena sativa (A. sativa) is a heterozygous hexaploid crop (AACCDD, 2n = 6x = 42) belonging to the genus Avena in the family Poaceae. It is extensively cultivated in 42 countries and territories worldwide due to its economic and nutritional significance [22], ranking as the sixth highest in global production [23]. A. sativa is recognized for its low carbon footprint and a wealth of health benefits, characterized by high concentrations of soluble fiber, β-glucan, lipids, proteins, and antioxidants [24]. Furthermore, A. sativa serves as a nutritious forage grass [25], demonstrates adaptability to diverse soil conditions, and shows significant drought tolerance [26], as well as considerable resistance to heavy metals [27]. Investigating the molecular mechanisms underlying drought tolerance and heavy metal resistance in A. sativa is essential for elucidating its stress tolerance strategies. The availability of genomic information for A. sativa enables in-depth exploration of the molecular mechanisms governing growth, development, and stress regulation at the genomic level [28, 29]. This research provides a foundation for a comprehensive study of the GST gene family and holds significant potential for guiding molecular breeding efforts aimed at enhancing resistance traits in A. sativa.

In this study, we identified the GST gene family members in A. sativa through bioinformatics analysis, focusing on chromosomal localization, phylogenetic relationships, conserved structural domains, gene structures, and gene duplication events. Additionally, we examined the expression profiles of four AsGST genes in roots subjected to PEG-induced drought stress and CdCl₂-induced cadmium stress at various time points, utilizing transcriptomic data and real-time quantitative PCR (RT-qPCR) technology. The results from both datasets exhibited a high degree of consistency, providing a robust basis for further validation of the functions of these four genes in response to drought and cadmium stress.

Results

Identification of GSTs in A. sativa

In this study, the protein sequences of 53 AtGST (Table S1) and 77 OsGST (Table S2) were utilized for the preliminary screening of the A. sativa genome via BLASTP and HMMER tools. Following further validation using CD-HIT software and the NCBI-CDD and SMART databases, a total of 118 non-redundant AsGST full-length genes, each containing complete GST_C and GST_N structural domains, were identified. The physicochemical properties of these AsGST proteins were assessed (Table S4), revealing protein lengths ranging from 178 amino acids (AsGSTF18) to 467 amino acids (AsEF1G4), with an average length of 248 amino acids. The molecular weights (MWs) of the AsGST proteins varied from 20.43 to 52.76 kDa, while isoelectric point (PI) values spanned from 4.79 to 7.82. Notably, five proteins (AsGSTU24, AsGSTF40, AsGSTF10, AsGSTF43, and AsGSTF3) exhibited PI values greater than 7, indicating a more basic nature compared to the other 173 proteins. Regarding subcellular localization, the majority of AsGST proteins were predicted to reside primarily in the cytoplasm and chloroplasts, with additional localization in the nucleus, extracellular matrix, peroxisomes, mitochondria, and vesicles. Furthermore, most AsGST proteins were classified as stable (instability index < 40); however, 46 proteins were identified as unstable, suggesting that while AsGST proteins are generally stable, some are prone to degradation. Additionally, only 21 proteins exhibited a hydrophobicity index greater than 0, indicating that most AsGST proteins are hydrophilic.

Phylogenetic analysis of GSTs in A. sativa

To investigate the evolutionary relationships among the GST proteins of AtGST, OsGST, and AsGST, a phylogenetic tree was constructed (Fig. 1). The AsGST gene family was classified into five distinct subfamilies based on the subfamily classifications of A. thaliana and Oryza sativa. Among these, the Phi subfamily contained the largest number of AtGST genes (59), followed by the Tau subfamily (50), while the Zeta and EF1Bγ subfamilies each had one member, as did the Lambda subfamily. Although subfamilies Theta, DHAR, and TCHQD are present in both A. thaliana and O. sativa, corresponding members were absent in A. sativa. Furthermore, no distinct clusters were identified for the Zeta and Lambda subfamilies, underscoring the considerable homology among AsGST members within each subfamily.

Fig. 1

Phylogenetic tree of GSTs from A. sativa (AsGST), A. thaliana (AtGST), and O. sativa (OsGST). These can be divided into eight subgroups (Tau, Lambda, Phi, Theta, TCHQD, Zeta, DHAR, and EF1Bγ); the different subgroups are distinguished by color

Analysis of the gene structure and motifs in A. sativa

To further investigate the interrelationships among the AsGST family genes, we predicted the conserved motifs of the AsGST proteins using the MEME suite, identifying a total of ten conserved motifs. Concurrently, the conserved structural domains of AsGST proteins were predicted with the NCBI-CDD tool. The distribution of the structural features, conserved motifs, and conserved structural domains of the AsGST genes is illustrated in Fig. 2. Results revealed that members of the same subfamily exhibit high similarity in conserved motif distribution, with motif 1 being the most conserved and, presenting across all members. With the exceptions of AsGSTL1 and AsGSTZ1, no significant differences in the number of conserved motifs were observed among the remaining GST sequences. Notably, specific motifs were exclusive to certain subfamilies; for instance, motif 7 was restricted to the Phi subfamily, while motifs 8 and 9 were unique to the Tau subfamily (Fig. 2B). These findings imply that members within the same subfamily may share functional similarities, whereas different subfamilies might display functional differentiation due to the specificity of their conserved motifs. A similar pattern emerged from the analysis of conserved structural domains, where GST sequences from the same subfamily possess highly similar domains, while unique structural domains were identified across different subfamilies. For example, only the EF1Bγ class GST genes contained the EF1G superfamily conserved structural domain (Fig. 2C), reinforcing the notion of functional similarity within subfamilies and potential functional differentiation between them. Additionally, we examined the coding sequences (CDS) and non-coding regions (UTRs) of the AsGST gene family (Fig. 2D). The number of exons varied from 1 to 11, with AsGSTL1 and AsGSTZ1 exhibiting the highest exon count of 11, while AsGSTU6, AsGSTU7, AsGSTU16, AsGSTU19, AsGSTU22, AsGSTU29, AsGSTU47, and AsGSTU49 each contained only one exon. Regarding introns, apart from AsGSTL1 and AsGSTZ1, which contained 8 introns, eight genes lacked introns entirely, and the remaining AsGST family members exhibited intron counts ranging from 1 to 7. Collectively, these results suggest notable differences in gene structure, conserved motifs, and structural domains among AsGST genes across different subfamilies, potentially contributing to variations in their physiological functions. Nonetheless, the overall genetic structure of the AsGST gene family appears conserved, particularly among members of the same subfamily, which may indicate shared functional roles.

Fig. 2

Phylogenetic analysis, conserved motifs, conserved domains, and exon/intron structure of AsGSTs. A A phylogenetic tree constructed based on the full-length sequences of AsGST proteins using the maximum likelihood method, with 1000 bootstrap replicates. B Conserved motifs 1–10 are represented by boxes of different colors. C Conserved domains are represented by boxes of different colors. D Non-coding regions are indicated by yellow boxes; exons are indicated by green boxes; introns are indicated by black lines. E Consensus sequence of each motif

Promoter cis‑acting element analyses in A. sativa

Using the PlantCARE online platform, we analyzed the cis-responsive elements present in the promoter regions of 118 AsGST gene sequences, extending approximately 2000 bp upstream of the transcriptional start site. The identified cis-acting elements in the AsGST gene promoters are depicted in Fig. 3. In total, 25 distinct cis-responsive elements were recognized within the AsGST gene family, classified into three functional categories: plant growth and development, abiotic and biotic stress responses, and phytohormone responses (Fig. 3B).

Fig. 3

Analysis of cis-acting elements in the promoter region of AsGST. A Visualization of cis-regulatory elements in the promoters of AsGST gene family members. B A heatmap displays the diversity of promoter elements in AsGST genes, with different colors and numbers representing different frequencies. C A histogram shows the cumulative counts of cis-acting elements in each category, represented by different colors

Within the plant growth and development category, G-box elements were the most abundant, constituting 51.03% of all identified elements and 17.49% of the total elements, primarily regulating transcription initiation frequency. Additionally, we identified A-box, CAT-box, O₂-site, and circadian elements linked to photorespiration and metabolic regulation; CCAAT-box and AT-rich elements associated with transcriptional regulation; RY-element and GCN4_motif elements related to seed development; as well as the HD-ZIP1 element pertinent to flower and fruit development. The diversity and complexity of these elements underscore the significant role of AsGST genes in plant growth and development. The second category encompasses elements related to abiotic and biotic stress responses, including ARE and GC-motif elements associated with anaerobic conditions, MBS elements tied to drought stress, LTR elements linked to low-temperature response, TC-rich repeats involved in plant defense mechanisms, and WUN-motif elements relevant to wound response. Notably, ARE elements were the most prevalent in this category, representing 36.50% of abiotic and biotic stress response elements and 6.33% of all elements, indicating the potential significance of AsGST genes in mediating plant responses to environmental stresses. The third category involves hormone-responsive elements, which include those for abscisic acid (ABRE elements), growth hormones (AuxRR-core and TGA-element elements), salicylic acid (CGTCA-motif and TCA-element elements), gibberellin (TATC-box, P-box, and GARE-motif elements), and jasmonic acid (TGACG-motif element). These hormonal elements play crucial roles in plant development and stress responses. Among these, ABRE elements were the most, accounting for 29.95% of hormone-responsive elements and 14.49% of all elements (Fig. 3C). Overall, the analysis of promoter regions highlights the multifunctionality of AsGST genes in plant growth and development, stress responses, and hormone signaling, operating through interactions with various cis-responsive elements.

Chromosome distribution and gene replication in A. sativa

Based on the annotation file of A. sativa, we mapped the chromosomal localization of AsGST genes (Fig. 4). The analysis revealed a heterogeneous distribution of the 118 AsGST genes across the 21 chromosomes, with a notable prevalence of genes located in both proximal and distal regions. Specifically, chromosome 4D (chr4D) harbored the greatest number of AsGST genes, with a total of 16, followed by chr1D with 12 genes and chr4 A with 11 genes. Conversely, chr6 C contained the fewest AsGST genes, with only 1 gene identified. This distribution pattern may be closely linked to the functional roles of AsGST genes. Gene duplication, through tandem or segmental mechanisms, is a common evolutionary phenomenon in plants that significantly contributes to genome amplification and diversification (Fig S1). To further investigate the evolutionary dynamics of the AsGST gene family, we analyzed the gene duplication events using the MCScanx tool (Fig. 5A). Our results indicated that 30 AsGST genes were involved in tandem duplication events, organized into 14 distinct clusters. Among these tandem duplications, 9 pairs were classified within the Phi subfamily, 4 pairs in the Tau subfamily, and 1 pair in the EF1Bγ subfamily. Additionally, we identified 109 segmental duplication pairs comprising 74 genes, suggesting that these duplication events may have played a crucial role in the expansion and functional diversification of the gene family. To assess the selective pressures acting on the duplicated genes, we calculated the ka/ks ratios. The results indicated that all gene pairs exhibited ka/ks values less than 1, ranging from 0.03 to 0.78, implying that purifying selection has been a significant force in the evolution of these AsGST gene pairs (Table S5).

Fig. 4

Distribution of AsGST on Avena sativa chromosomes. The chromosomal positions of each AsGST gene were mapped onto the A. sativa genome

Fig. 5

Collinearity analysis of AsGST. A Intraspecific collinearity analysis of AsGST (paralogous gene pairs are highlighted by red lines. Gene density distribution on chromosomes is represented by a heatmap and dashed lines.) B Collinearity analysis of AsGST in Arabidopsis thaliana, Oryza sativa, Triticum aestivum, and Avena sativa (collinear GST gene pairs are indicated by red lines, while collinear regions between genomes are indicated by gray shading)

To further elucidate the phylogenetic mechanisms and conserved genomic structure of the AsGST gene, we examined the covariance between A. sativa with A. thaliana, O. sativa, and Triticum aestivum (Fig. 5B). Our analysis revealed a single pair of segmental duplications in the GST gene between A. sativa and A. thaliana, whereas there were 36 and 63 pairs of segmental duplications identified between O. sativa and T. aestivum, respectively. This disparity highlights a closer genomic relationship between A. sativa and both O. sativa and T. aestivum, particularly in the expansion and conservation of the GST gene family. These findings suggest shared molecular mechanisms that facilitate adaptation to environmental stresses and evolutionary processes among these species. The observed conservatism may stem from the common ancestry and similar genomic evolutionary histories shared by A. sativa, O. sativa, and T. aestivum. In contrast, A. thaliana, as a dicotyledonous crucifer, exhibits a relatively distant genetic relationship with A. sativa.

Expression of AsGST under osmotic/drought stress conditions in A. sativa

In our investigation into the involvement of AsGST genes in the drought stress response of A. sativa, we analyzed changes in the expression of 118 AsGST genes in the root system of A. sativa using transcriptomic data sourced from the NCBI database (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA1056521/). Our analysis revealed a lack of detectable expression in 13 of these genes across all treatment times, while 83 genes showed no significant variation in expression (Fig. 6A, Table S6). Notably, 22 genes exhibited significant alterations in expression, among them 15 from the Tau subfamily (AsGSTU6 to AsGSTU18, AsGSTU21, AsGSTU37, AsGSTU43, AsGSTU45) and seven from the Phi subfamily (AsGSTF12, AsGSTF21, AsGSTF24, AsGSTF33, AsGSTF36, AsGSTF57, AsGSTF58). In particular, AsGSTU37 showed marked down-regulation for drought (PEG) treatment (osmotic stress mimicking drought conditions in the soil), especially at the initial stages (6 h and 12 h) where expression levels dropped by approximately threefold and sixfold. The other 21 genes showed significantly up-regulated expression (Table S6). Following these findings, we selected four genes showing similar expression patterns for confirmation via real-time quantitative PCR (RT-qPCR) from samples exposed to PEG treatment. To bolster the reliability of our experiment, in addition to validating the root transcriptome data, we quantified the gene expression levels in the leaves of A. sativa, and extended the treatment duration to 48 h. The RT-qPCR results (Fig. 7B) revealed dynamic, generally increasing expression patterns in roots for AsGSTU12, peaking at 12 h and 48 h with about tenfold up-regulation, though at 24 h, the expression fell but remained approximately threefold above control levels. Similarly, AsGSTU13, AsGSTU14, and AsGSTU15 showed trends of initial increase followed by a decrease, with AsGSTU13 expression levels rising from around twofold up-regulation at 6 h to 16-fold at 12 h, before gradually falling to fivefold at 72 h (Fig. 7B). Conversely, in A. sativa leaves, the expression patterns of these four genes differed (Fig. 7A). Significant up-regulation was observed in the early stages of drought treatment (6 h or 12 h), but expression levels significantly declined by the mid-point of the drought treatment (12 h), to levels on par with, or even lower than, the untreated condition (Fig. 7A). After being treated with PEG for 72 h, AsGSTU12, AsGSTU13, and AsGSTU15 were up-regulated about twofold, and AsGSTU14 was up-regulated about fivefold (Fig. 7A).

Fig. 6

Transcript abundance map of AsGST. A A heatmap shows the expression profiles of AsGST genes in roots under drought stress; B A heatmap shows the expression profiles of 22 differentially expressed genes in roots under drought stress; C A heatmap shows the expression profiles of AsGST genes in leaves and roots under cadmium stress; D A heatmap shows the expression profiles of 4 differentially expressed genes in leaves under cadmium stress; E A heatmap shows the expression profiles of 11 differentially expressed genes in roots under cadmium stress

Fig. 7

RT-qPCR analysis of four differentially expressed genes. A Relative expression levels of four genes in leaves after PEG treatment; B Relative expression levels of four genes in roots after PEG treatment; (C) Relative expression levels of four genes in leaves after CdCl₂ treatment; D Relative expression levels of four genes in roots after CdCl₂ treatment. The Y-axis and X-axis represent the six time points of treatment and relative expression levels, respectively. The relative expression levels of the genes at 0 h were set to 1 and calculated using a normalization method. The mean ± standard deviation (SD) was obtained from three biological replicates and three technical replicates. Error bars indicate the standard deviation. Different letters indicate significant differences at P < 0.05 level, while the same letters indicate no significant differences

Expression of AsGST under cadmium stress in A. sativa

To assess the function of AsGST genes in response to cadmium stress, we examined transcript abundance alterations of 118 AsGST genes in A. sativa leaves and roots under cadmium stress using transcriptome data (https://www.ncbi.nlm.nih.gov/bioproject/ PRJNA1116317) (Fig. 6C). In both root and leaf, 13 AsGST genes remained undetected. A total of 101 AsGST genes in leaves and 94 in roots displayed non-significant differential expression patterns as indicated in Table S7. As shown in Fig. 6D and 6E, four genes demonstrated significant expression differences in leaves and 11 in roots, respectively. Notably, AsGSTU12, AsGSTU13, AsGSTU14, and AsGSTU15 displayed significant upregulation (4–13 folds) in both leaves and roots in response to cadmium treatment (Fig. 6D-E). To validate these results, we selected genes significantly upregulated in both roots and leaves (AsGSTU12, AsGSTU13, AsGSTU14, AsGSTU15), and assessed their expression using RT-qPCR across four time points (6 h, 12 h, 48 h, and 72 h). RT-qPCR results showed early up-regulation (within 12 h) in leaves of all four genes in response to cadmium treatment, but no significant difference or even suppression was detected afterward time points (Fig. 7C). In A. sativa roots, a transient up-regulation of four genes was observed, reaching a peak at 6 h of cadmium treatment, and then declining, with a 30-fold increase for AsGSTU12, 27-fold for AsGSTU13, 45-fold for AsGSTU14, and 20-fold for AsGSTU15 (Fig. 7D). Although the RT-qPCR results showed varying degrees of up-regulation compared to the transcriptome data, the overall trend remained consistent.

Subcellular localization analysis of AsGSTU15 protein

To elucidate the expression pattern of GST proteins in A. sativa, the AsGSTU15 protein, with a significant response to both drought and cadmium stress, was chosen for the subcellular localization analysis. The AsGSTU15 fused with the green fluorescent protein (35S::AsGSTU15-GFP) was transformed in tobacco leaf, and expressed in both cytoplasm and nucleus of tobacco epidermal cells, suggesting the localization of AsGSTU15 protein in these regions (Fig. 8). This finding aligns with the prediction of the Wolf website, which postulated the cytoplasmic localization of AsGSTU15 protein primarily (Table S4).

Fig. 8

Subcellular localization of AsGSTU15. Bars = 20 μm

Upstream regulatory mechanism of AsGSTU gene and prediction of interacting proteins

To explore the regulatory mechanisms of AsGSTU12, AsGSTU13, AsGSTU14, and AsGSTU15 genes in response to drought and cadmium stress, we postulated their potential upstream transcription factors (TFs) and interacting proteins. Utilizing transcriptome data from A. sativa roots under these stress conditions, we identified 394 common DEGs (Fig. 9A). From this differential gene set, a potential transcription factor regulatory network was conceived (Fig. 9B). Our analysis revealed six transcription factors that may potentially regulate these target genes. Notably, two genes from the ERF family, AVESA.00010b.r2.2 CG0297950 and AVESA.00010b.r2.2DG0331590, demonstrated a high likelihood of binding to AsGSTU15 and AsGSTU12, with their expression levels notably amplified under drought and cadmium stresses (Fig. 9D and E). These genes may act as key upstream transcription factors. Additionally, AVESA.00010b.r2.3DG0552900 and AVESA.00010b.r2.4 AG0632530 from the NAC transcription factor family, which also demonstrated elevated expression levels under the stress conditions, may exert some regulatory influence. The prediction of interacting proteins revealed that from the differential gene set, AsGSTF58, from the same GST gene family, had the closest interaction with the target genes and its expression level was amplified under stress conditions. This was followed by AVESA.00010b.r2.1DG0134140, which is homologous to AtBGLU, and AVESA.00010b.r2.6 CG1099650 gene, homologous to AtGER5, both identified as potential interacting proteins for future examination.

Fig. 9

Prediction of upstream transcription factors and interacting proteins for four differentially expressed genes, and the expression levels of different genes in the transcriptome under drought and cadmium stress. A Venn diagram of different genes in roots under cadmium and drought stress. The two gene sets are filled with different colors, and the overlapping area represents the genes shared between the gene sets, with numbers indicating the count of shared genes. B, C Prediction of upstream transcription factors and interacting proteins for the target genes. Nodes represent gene or protein names, with pink nodes indicating transcription factors. The size of the nodes represents gene connectivity, with larger nodes indicating stronger connectivity. Lines connecting nodes represent the interaction relationships between genes or proteins. The color of the lines (min p-value) represents the probability of motif occurrence or protein interaction, with smaller p-values indicating more reliable results. The thickness of the lines (max score) represents the motif occurrence score, with thicker lines indicating a higher possibility of transcription factor binding to the output sequence or stronger protein interaction. D A heatmap shows the expression profiles of six predicted transcription factors in roots under drought stress. E A heatmap shows the expression profiles of six predicted transcription factors in roots under cadmium stress. F A heatmap shows the expression profiles of four predicted interacting proteins in roots under drought stress. G A heatmap shows the expression profiles of four predicted interacting proteins in roots under cadmium stress

Physiological response of A. sativa to drought and cadmium stress

To investigate the physiological responses of A. sativa to drought and cadmium stress, we analyzed dynamic changes in GSH content in shoots and roots under different stress durations (Fig. 10). Under drought stress (Fig. 10B, C), shoot GSH content increased significantly over time, reaching the highest at 48 h. In contrast, root GSH content decreased after an initial reduction at 6 h, with the lowest value at 72 h. Under cadmium stress, shoot GSH demonstrated an initial elevation at 6 h, followed by successive fluctuations characterized by a secondary peak at 48 h before declining to baseline levels at 72 h (Fig. 10D). This oscillatory response differed markedly from the consistent depletion observed in cadmium-stressed roots (Fig. 10E), where glutathione content exhibited progressive reduction from 6 h onward, maintaining persistently low concentrations through 48–72 h of treatment.

Fig. 10

Growth status and GSH content of Avena sativa. A The growth status of A. sativa at six time points (0–72 h); the scale bar is 10 cm. B The analysis results of GSH content. The mean ± standard deviation was obtained from six biological replicates. Error bars indicate the standard deviation. Different letters indicate significant differences, while the same letters indicate no significant differences at P < 0.05 level

Discussion

The GST gene family, characterized by its multifunctionality, is ubiquitous in plant species and has gathered considerable interest due to its functional complexity. Comprehensive genome-wide analyses have been conducted across an array of plants, including A. thaliana [30], O. sativa [31], S. lycopersicum [32], T. aestivum [33], Hordeum vulgare [34], G. hirsutum [13], Glycine max [35], and Solanum tuberosum [36]. Notably, the sequencing of the A. sativa genome has been successfully accomplished [28, 29], promoting advanced investigations into the function of the GST gene family.

Identification and bioinformatics analysis of the AsGST gene family

In our research, we conducted a thorough characterization of the GST gene family in A. sativa. The bioinformatic analysis identified 118 AsGST genes, classified into five distinct subfamilies: Tau, Phi, Zeta, Lambda, and EF1Bγ. This classification was based on the subfamily categorization in A. thaliana and O. sativa. Phylogenetic examinations revealed that the AsGST members closely resembled those of the same GST category in A. thaliana and O. sativa—a similarity that suggests a shared ancestral lineage. Interestingly, despite the considerable size of the A. sativa genome (10.76 Gb) [29], which is approximately 79.6 times larger than that of A. thaliana (135 Mb) [37] and about 40 times that of O. sativa (385.7 Mb) [38], the AsGST genes did not exhibit the Theta, DHAR, and TCHQD subfamilies observed in A. thaliana and O. sativa. This indicates that the AsGST's evolutionary trajectory may involve gene loss or divergence, associated with functional substitution and environmental adaptation [39], and may not be directly influenced by genome size. Of all the identified classes, the Phi subfamily, with 59 members, was the largest in A. sativa, followed by the Tau subfamily with 40 members. This finding aligns with the widespread and distinct distribution of Tau and Phi subfamilies in plant species [40]. These genes have shown importance in plant response to various stresses, including salt [41], oxidative stress [42], heavy metal, and drought [43]. During species evolution, gene duplication is essential for the development of new biological functions and the expansion of gene families [42]. Gene family expansion is mainly realized by both segmental and tandem duplication [43]. Our research identified 108 segmental duplication pairs and 14 tandem duplication clusters unevenly distributed across the 21 chromosomes of the AsGST family, which are crucial for the evolution of the AsGST gene family. In this study, segmental duplication phenomena were predominantly observed, secondary to the role of tandem duplications, a trend similarly noted in the GST family of T. aestivum [43]. Moreover, the ka/ks ratios for all duplicated AsGST gene pairs were less than 1, suggesting that comprehensive purifying selection has occurred in all AsGST genes. This mechanism is likely critical for preserving gene functionality and promoting adaptability to environmental stressors.

The configuration of exon–intron structures is integral to gene evolution. Upon scrutinizing gene structures and conserved motifs, it was observed that varying GST types possess distinct motifs and gene structures. However, AsGST genes from the same subfamily display similarities in terms of exons and introns count and positioning, and they share comparable motifs. This finding implies potential expansive functional redundancy within the AsGST family, and suggests that AsGST genes from an identical subfamily may perform similar biological functions in A. sativa. Moreover, the number of exons/introns within AsGST family members ranged from zero to eleven, possibly indicative of deletion or insertion events throughout the evolutionary history of A. sativa [44].

Members of the GST family can respond to a broad spectrum of abiotic and biotic stresses. Regulatory elements known as cis-acting elements are integral to plant adaptation to these stresses, alongside growth and development, by controlling the expression of pertinent genes [45]. Analysis through PlantCARE revealed that the AsGST gene's promoter region encompasses three categories of cis-acting elements: those responsive to plant growth and development, abiotic and biotic stress, and phytohormone. These cis-acting elements may influence the expression level and response capability of AsGST family genes under various conditions such as growth, development, stress response, and hormone regulation. The results from the cis-acting element analysis suggest the potential involvement of AsGST genes in the regulation of A. sativa growth and development processes, such as seed, flower, and fruit development, and in responses to various abiotic and biotic stresses, which could include drought and heavy metal exposure. They may also be subject to regulation by hormones, such as abscisic acid and growth hormone, among others. The diverse array of cis-acting elements in the promoters of AsGST suggests potential synergistic interactions among the promoters, collectively enhancing the responsiveness of AsGST to a range of stimuli. This intricate regulatory network might enable AsGST to adapt flexibly to shifting environmental conditions, thus conferring stress (eg. drought and heavy metal) tolerance to A. sativa.

The predicted upstream TFs with the highest binding probability were AVESA.00010b.r2.2 CG0297950 and AVESA.00010b.r2.2DG0331590, both belonging to the ERF family and functioning as ethylene response factors in plants. BLASTX analysis revealed that AVESA.00010b.r2.2 CG0297950 is analogous to the AtERF114, which plays a crucial role in plant regeneration and development by regulating growth hormone signaling in response to mechanical and hormonal signals, thereby influencing the developmental and regenerative organ responses [46]. AtERF114 also controls the synthesis of secondary metabolites by upregulating AtPAL1 transcription, consequently enhancing resistance to fungal pathogens in A. thaliana [47]. AVESA.00010b.r2.2DG0331590, conversely, is analogous to AtERF113, which plays a significant role in plant response to abiotic stresses. For instance, GmERF113 was found to positively regulate drought response by activating GmPR10 - 1 expression, thereby improving plant drought resistance [48]. Moreover, RhERF113, induced by ethylene in floral organs and upregulated during flower senescence, can regulate cytokinin content during ethylene-induced petal senescence [49]. In A. sativa, these two genes may enhance plant tolerance to drought and cadmium stresses by regulating AsGSTU12 or AsGSTU15 expressions. Subcellular localization experiments indicated that the AsGSTU15 protein is located in both the cytoplasm and nucleus (Fig. 8), and given that transcription factors ordinarily function in the nucleus, ERF family transcription factors AVESA.00010b.r2.2 CG0297950 and AVESA.00010b.r2.2DG0331590 may regulate AsGSTU15 expression in the nucleus. Protein interaction predictions showed that AsGSTF58 has the highest interaction probability with AsGSTU12, AsGSTU13, AsGSTU14, and AsGSTU15 (Fig. 9). Although these proteins belong to different GST subfamilies, their structural and functional similarities suggest joint participation in antioxidant response and heavy metal detoxification processes in response to drought and cadmium stresses. Subsequently, the proteins with the second highest binding probability were AVESA.00010b.r2.1DG0134140 and AVESA.00010b.r2.6 CG1099650. BLASTX comparison showed that AVESA.00010b.r2.1DG0134140 is analogous to the AtBGLU6, a β-glucosidase gene of the glycoside hydrolase family 1 (GH1) involved in various biochemical processes [50]. In O. sativa, the BGLU6 gene (Os3BGlu6) can hydrolyze glucose-bound ABA (ABA-GE) to free ABA, promoting stomatal closure, reducing water evaporation, and enhancing plant drought tolerance [51]. Additionally, AVESA.00010b.r2.6 CG1099650 is homologous to the AtGER5 gene, another ABA and stress response protein involved in ABA-mediated stress response. Increased ABA levels under abiotic stress conditions, such as drought, correspond to an upregulation of AtGER5 expression, aiding plants in regulating their response to drought stress [52]. Hence, AtBGLU6 and AtGER5 homologs in A. sativa may functionally complement these four AsGSTU genes and synergistically participate in signaling and regulation, thereby enhancing plant tolerance to abiotic stress and exhibiting higher binding probability with target genes in the protein interaction network.

The AsGSTU subfamily may respond to PEG-induced drought stress by regulating redox homeostasis through scavenging ROS, participating in the ABA signaling pathway, and regulating osmotic homeostasis

Drought stress substantially impedes plant growth and crop yield, with the expression pattern of the GST gene family closely associated with plant resistance [32]. GSH is a critical regulator of cell growth and proliferation, with alterations in its content closely linked to a plant's physiological state [53]. The overexpression of AtGST17 in A. thaliana led to a marked increase in the GSH/GSSG ratio and a greater accumulation of GSH and ABA, thereby regulating the intracellular redox balance, reducing stomatal aperture, lowering transpiration rates, and ultimately enhancing drought tolerance [54]. The primary sites for ROS production are chloroplasts and mitochondria within above-ground leaves [55], and GSH plays a crucial role in mitigating ROS damage [56]. For instance, overexpression of PtrGSTU23 in Populus trichocarpa [57] and CsGSTU8 in A. thaliana [18] reduced in vivo ROS content by catalyzing reactions between GSH and ROS, thereby improving drought tolerance. Under drought stress, in vivo ROS such as H₂O₂ content in A. sativa significantly accumulated, leading to a consequential rise in GST enzyme activity [58]. The present study found that A. sativa growth and development were markedly hindered under drought conditions (Fig. 10A), possibly due to the altered resource allocation under stress conditions, where resources were primarily used to manage drought stress instead of growth and development. Changes in A. sativa's GSH content were intimately associated with its physiological status. Under drought stress, shoot GSH content in A. sativa significantly decreased in the initial stage (6 h), potentially due to the sharp increase in ROS levels at this time, using large amounts of GSH for ROS scavenging to protect cells from oxidative stress damage. Over time, the GSH content gradually increased, peaking at 48 h, suggesting that A. sativa may have initiated corresponding regulatory mechanisms, such as improving the antioxidant defense system and increasing the synthesis or recycling of GSH to maintain intracellular GSH homeostasis and more effectively manage persistent ROS stress. Despite a decrease at 72 h, the GSH content remained significantly higher than the control, indicating A. sativa's strategy to maintain redox homeostasis under prolonged drought stress (Fig. 10B). Conversely, root GSH content progressively declined under drought stress, reaching a minimum at 72 h (Fig. 10C). This may suggest that roots experienced substantial inhibition under drought stress (Fig. 10A) and, as they are not the primary ROS-producing organ, may not require excessive GSH for ROS scavenging.

Transcriptome data analysis highlighted that 22 AsGST genes were markedly differentially expressed in A. sativa roots under drought treatment, most of which belonged to the Tau subfamily. Notably, four genes (AsGSTU12, AsGSTU13, AsGSTU14, and AsGSTU15) were significantly upregulated, with their transcript levels increasing approximately threefold under drought stress. The RT-qPCR validation results, consistent with the transcriptome data, affirmed the positive role of these four genes in the drought stress response. BLASTX comparison revealed that AsGSTU12 shared the highest homology with AtGSTU15, AsGSTU13 and AsGSTU15 with AtGSTU17, and AsGSTU14 with AtGSTU18. No direct studies have elucidated the specific functions of AtGSTU15 and AtGSTU18 yet. The prediction of promoter cis-acting elements suggests that all four AsGSTU genes of A. sativa may contain ABRE elements (abscisic acid responsive elements), and it has been determined that AtGSTU17 expression is induced by ABA, which negatively regulates plant drought tolerance [54]. Under drought conditions, atgstu17 displayed a lower water loss rate and smaller stomatal opening compared to the wild type, demonstrating greater drought tolerance. Concurrently, atgstu17 accumulated higher levels of GSH under normal growth conditions, aiding in maintaining cellular redox balance under adverse conditions [54]. Additionally, AtGSTU17 is implicated in the expression of hypocotyl elongation, anthocyanin accumulation, and far-red light-mediated greening inhibition in A. thaliana seedlings, with its expression level regulated by various photoreceptors, especially photoreceptor phytochrome A (phyA), which functions under all light conditions [59]. Therefore, it is hypothesized that the AsGSTU13 and AsGSTU15 genes (homologous of AtGSTU17) may exhibit similar functions in response to drought stress in A. sativa. These genes may be both regulated by ABA signaling and express fluctuations in up- or down-regulation in A. sativa leaves and roots under different temporal drought treatments (Fig. 7A and B) controlling GSH content (Fig. 10B, C) for scavenging excess ROS to avert oxidative damage and ensure a basic level of ROS signaling function [58]. However, AsGSTU12, AsGSTU13, AsGSTU14, and AsGSTU15 were upregulated in roots at various treatment times because roots are the first to sense drought and water deficit, and persist in these conditions. It is posited that the GSTU gene family effectively promotes the GSH-catalyzed detoxification to scavenge ROS and prevent root over-oxidization, hence the GSH content decreased (Fig. 10C). The AtGSTU18, homologous to AsGSTU14, has been shown to scavenge the reactive carbonyl compound acrolein (RCS) produced by ROS oxidizing the plasma membrane, thus mitigating the toxic effects of RCS on cells [60]. There are no known studies on the AtGSTU15, homologous to the AsGSTU12. However, its significant upregulated expression in A. sativa roots and transient expression in leaves may be similarly induced or regulated by ABA, as well as involved in ROS scavenging and signaling roles.

Osmoregulatory substances like proline and betaine are integral for preserving osmotic equilibrium and safeguarding cellular structure and function [61]. Notably, contents of these substances significantly escalated in Nicotiana tabacum with overexpressed GSTU gene under drought stress conditions [62]. Conversely, osgstu17 displayed a significant reduction in leaf proline content, highlighting the crucial role of OsGSTU17 in positively regulating proline synthesis [63]. Consequently, we conjecture that the four AsGSTU genes in A. sativa may contribute to osmotic regulation under PEG-induced drought stress, thereby maintaining intracellular osmotic pressure equilibrium and cellular membrane stability, and ultimately augmenting the plant's drought resilience.

The AsGSTU subfamily may respond to cadmium stress by regulating redox homeostasis through scavenging of ROS and promoting cadmium chelation and compartmentalization

Cadmium, an extremely toxic heavy metal pollutant, significantly disrupts plant growth and development. Unlike drought stress, cadmium exerts direct toxicity on plant cells, severely hampering A. sativa growth and development following cadmium treatment (Fig. 10A). Cadmium stress results in an excessive ROS accumulation, thus inducing oxidative stress. This process involves production and signal transduction of various substances, which cause oxidative damage to biomolecules (eg. proteins, lipids, and nucleic acids), subsequently affecting enzymatic activity and cellular membrane integrity [64]. The GST gene family is pivotally involved in plant response to cadmium stress. The GST enzyme, by facilitating the combination of GSH with the cadmium ion to form either non-toxic or less toxic complexes, diminishes cellular cadmium toxicity, promotes its detoxification and sequestration, and thus enhances plant tolerance to cadmium stress [65].

Under cadmium stress, the GSH content in the shoot of A. sativa exhibited a significant fluctuation, reflecting its role in ROS scavenging and cadmium ion chelation (Fig. 10D). At early stage of the cadmium stress (6 h), a substantial increase in GSH content occurred alongside the upregulation of AsGSTU12, AsGSTU13, and AsGSTU15 expressions (Fig. 7C). GSH, functioning as a crucial antioxidant, rapidly scavenges excess ROS thereby mitigating cellular oxidative damage. Additionally, GSH forms stable complexes with cadmium ions, effectively immobilizing them within cells and/or transporting them to vesicles, thus reducing cytoplasmic cadmium ion concentration and preventing organellar damage caused by cadmium. However, as cadmium stress continues, GSH content fluctuates and corresponding AsGSTU14 and AsGSTU15 show differential upregulation (Fig. 7C). This response likely reflects a balance of ongoing ROS scavenging, maintenance of a certain ROS level for signaling purposes, and compartmentalization of cadmium ions facilitating sustained, albeit slow, growth and development (Fig. 10A). Similarly, Cosmos bipinnatus enhances cadmium tolerance via the augmentation of antioxidant levels attributed to an increase in GSH content [66]. O. sativa strains that overexpress the OsGST (LOC_Os10 g38160) gene display augmented cadmium resilience and diminished cadmium accretion, inferring that the OsGST gene modulates cadmium accumulation and resistance in rice through the preservation of redox equilibrium [67]. It has also been observed that overexpression of either OsGSTU5 or OsGSTU37 lines significantly elevated GSH content and enhanced seed germination, seedling growth, and survival under cadmium-related stress [65]. Consistent with this, overexpressing of OsGSTU6 resulted in a decreased cadmium accumulation in foliage, thereby boosting plant tolerance to cadmium-induced stress [68]. Conversely, under cadmium stress, the GSH content in A. sativa roots decreases, reaching a nadir at 48–72 h (Fig. 10E). This decrease likely reflects cadmium-induced metabolic disruption in roots, where GSH is prioritized for ROS scavenging and cadmium chelation. Despite the depletion of GSH, GSTU family genes (AsGSTU12, AsGSTU13, AsGSTU14, and AsGSTU15) were upregulated after 6 h of cadmium treatment, with most exhibiting significantly higher expression levels than controls at subsequent time points. This suggests that the heavy consumption of GSH for ROS scavenging and continual cadmium ion chelation depletes GSH content, and potentially inhibits the activity and gene expression of enzymes crucial to GSH synthesis and metabolism. Overexpression of PuGSTU17 enhances resistance to Zn²⁺ stress by increasing root GST content and reducing H₂O₂ content [69]. The expression profile of GSTU family genes in different parts of the plant over time reinforces their role in ROS scavenging and cadmium ion binding, underscoring their critical function in A. sativa's resistance to heavy metal cadmium.

Conclusions

In this study, an in-depth investigation was conducted on the composition, structure, and stress responses of the AsGST gene family. Bioinformatic analysis identified 118 AsGST genes grouped into five subfamilies, which exhibited significant fluctuation in their expression patterns under adverse stress. Notably, four genes within the Tau subfamily (AsGSTU12, AsGSTU13, AsGSTU14, and AsGSTU15) demonstrated significant upregulation under both drought and cadmium stress. It was postulated that these four AsGSTU genes contribute towards the catalytic scavenging of ROS via GSH substrates, thus regulating redox equilibrium and maintaining basic ROS signaling roles. Under cadmium stress, these genes are hypothesized to respond by facilitating the scavenging of ROS using GSH to regulate redox homeostasis, promoting the chelation and fixation of GSH with cadmium, or by compartmentalizing cadmium under heavy metal stress. This study revealed the conjoint effects of these four AsGSTU genes under drought and cadmium stress, offering valuable genetic resources for A. sativa stress tolerance breeding and adversity adaptation research.

Materials and methods

Identification and classification of GST genes from A. sativa

In this research, the genome and annotation data of A. sativa were procured from the GrainGenes database (https://wheat.pw.usda.gov/GG3/content/avena-sang-download). Concurrently, the protein sequences of 53 AtGST (Table S1) were retrieved from the TAIR database (https://www.arabidopsis.org/) [30]. From the Ensembl plant database (http://plants.ensembl.org/index.html), we obtained the protein sequences of 77 OsGST (Table S2), although the protein sequences for OsGSTU3 (LOC_Os10 g38501) and OsGSTU4 (LOC_Os10 g38495) were not found [31]. These 130 assembled sequences were utilized as reference protein sequences for a BLASTP search employing TBtools-II software (version 2.119) [70], with an E-value ≤ 1e⁻¹⁰ and identity ≥ 50% as screening criteria, to identify GST genes in A. sativa [71]. The subsequent screening of the potential AsGST genes was performed using HMMER software (version 3.0), applying the Hidden Markov Model (HMM) of characteristics of GST-N (PF02798) and GST-C (PF00043) structural domains from the PFAM database (https://pfam-legacy.xfam.org/), with an E-value threshold set at < 1 × 10^⁻⁵. Redundant sequences from these protein sequences were removed using CD-HIT (version 4.8.1) with parameters set at c = 0.9 and n = 5. Only the longest transcripts were selected for AsGST genes with multiple transcripts [72]. The integrity of conserved structural domains in the protein sequences of candidate genes was examined using the NCBI conserved domain database (https://www.ncbi.nlm.nih.gov/cdd/) and the Simple Modular Architecture Research Tool (SMART) website (https://smart.embl.de/), discarding those potential genes not containing GST-conserved structural domains. Physicochemical properties of AsGST proteins, including protein length, molecular weight, isoelectric point, instability index, and average hydrophobicity, were analyzed via the ExPASy-ProtParam website (https://www.expasy.org/resources/protparam) [73]. Predictions of the subcellular localization of the AsGST were made using the WOLF PSORT website (https://wolfpsort.hgc.jp/).

Chromosome localization and phylogenetic tree construction

Utilizing TBtools-II software gene localization visualization tool, we ascertained the chromosomal locations of different AsGST gene classes using annotation data from A. sativa [70]. To delve into the phylogenetic relationships and evolutionary history of AsGST genes, we amalgamated the identified AsGST proteins with 150 GST proteins from A. thaliana and O. sativa referenced in prior studies. This was followed by a multiple sequence comparison using the default parameters of Clustal X software. Based on the categorization of A. thaliana and O. sativa subfamilies in the sequence comparison results and the structural features of AsGST proteins, we classified and named the AsGST genes. In A. sativa, gene names commence with the prefix “As” followed by a category identifier. Hence, the AsGST genes were ultimately sorted into five subfamilies: AsGSTU, AsGSTF, AsGSTZ, AsGSTL, and AsEF1BG, representing the Tau, Phi, Zeta, Lambda, and EF1Bγ subfamilies, respectively [74, 75]. Lastly, we constructed a neighbor-joining phylogenetic tree encompassing all these genes using MEGA software (version 11.0) and employed the Interactive Tree of Life (iTOL) website (https://itol.embl.de/) for phylogenetic tree refinement.

Analysis of the conserved domain and gene structure

We conducted a conserved motif analysis of AsGST proteins using the MEME online platform (https://meme-suite.org/meme/). In this analysis, we aimed to identify 10 conserved motifs, keeping other parameters at default settings. Following this, we inputted the phylogenetic tree file, conserved motif information, genome annotation data, and the protein structure domain prediction file derived from the NCBI-CDD tool (https://www.ncbi.nlm.nih.gov/cdd/) into the gene structure view component of the TBtools-II software. This comprehensive procedure facilitated the visualization and analysis of the phylogenetic tree, conserved motifs, structural domains, and gene structures of AsGST genes.

Cis‑acting element analysis of AsGST promoter region

Utilizing the annotation file of A. sativa, sequences 2000 bp upstream of the start codon of the AsGST gene were procured through the GTF/GFF3 sequence extraction instrument of TBtools-II. These sequences underwent analysis using the PlantCARE database (https://bioinformatics.psb.ugent.be/webtools/plantcare/html/) to discern cis-acting elements within the promoter region of the AsGST gene. Subsequently, the promoter region of the AsGST gene underwent visualization and enumeration of the cis-acting elements within its promoter region by use of the Simple BioSequence Viewer tool of TBtools-II. The fresh data were then collated and displayed in Excel for further scrutiny.

Analysis of gene duplication and collinearity

We employed the One Step MCScanX tool within TBtools-II software, setting the E-value limit below 1 × 10^⁻¹⁰, to identify AsGST gene duplication events. The Advance Circos program was subsequently used to visualize and analyze the location, gene density, and covariance details of the AsGST gene. Concurrently, we employed multiple covariate scanning methods to further identify homologous genetic linkages within the A. sativa genome and across other species. Additionally, the Simple Ka/Ks Calculator tool in TBtools-II software was utilized to calculate the ratio of non-synonymous to synonymous mutation rates (Ka/Ks) for two protein-coding genes, which is a pivotal determinant for the selection pressure on the genes.

Plant material and stress treatment

A. sativa seeds were initially screened using a crop seed tester (BIO-seed M-P, Research Center of Information Technology, Beijing Academy of Agriculture and Forestry Sciences, Beijing, China) based on the length, width, and perimeter of the seeds. The seeds then underwent sterilization via a 20-min immersion in a 1.5% sodium hypochlorite solution, prior to being relocated to germination boxes. A. sativa seedlings were grown in a light-shielded growth chamber at a stable temperature of 25 ± 1 °C until attaining a height between 2–3 cm. Seedlings were then transferred to the two-true leaf phase in a 1/2 Hoagland nutrient solution that contained either 15% (w/v) polyethylene glycol 6000 (PEG 6000) or 100 μM CdCl₂ for subsequent treatment. Following different treatment durations at 0 h, 6 h, 12 h, 24 h, 48 h, and 72 h, samples were collected. The shoot and root portions of plants were cleansed with sterile water, excess water was eliminated, and they were quickly frozen using liquid nitrogen. The samples were then preserved at − 80 °C in a refrigerator for future experimental analysis.

RNA extraction and RT‑qPCR

Total RNA extraction from the leaves and roots of A. sativa was conducted at each time point utilizing an RNAiso Plus kit (Takara, Japan) as per the kit's user manual. Subsequently, Nanodrop spectrometry (Thermo Fisher Scientific, USA) was employed to ascertain the concentration and purity of the total RNA. The PrimeScript™ RT kit (Takara, Japan) was then used to synthesize first-strand complementary DNA (cDNA), using 1 µg of total RNA as a template. Furthermore, primer design for RT-qPCR experiments (Table S3) was executed using Primer Premier software (version 5.0), and primer specificity was confirmed using the Primer Check function of TBtools-II software. RT-qPCR assays were carried out using the CFX Connect instrument (Bio-Rad, USA) and ArtiCan^CEO SYBR qPCR Mix reagents (Tsingke, China). The total reaction system was set to a volume of 20 µL, comprising 10 µL of ArtiCan^CEO SYBR qPCR Mix, 0.4 µL each of upstream and downstream primers (at 10 µM concentration), and 2 µL of template cDNA, diluted tenfold with water. The remaining volume was topped up to 20 µL with ddH₂O. The RT-qPCR reaction parameters were as follows: Pre-denaturation at 95 °C for 5 min, followed by denaturation at 95 °C for 10 s, and annealing at 60 °C for 30 s, with a total of 40 cycles. The default values of the instrument were used for the melting curve acquisition program. The internal reference gene GAPDH was used to normalize cDNA concentration [76]. Lastly, the relative gene expression levels were determined using the 2^−ΔΔCT method [77].

Subcellular localization

We integrated the CDS of the AsGSTU15 gene, which lacks a stop codon, into the 35S::GFP binary vector based on pBWA(V)HS via molecular cloning techniques. Then we electroporated the constructed plasmid into the GV3101 strain of Agrobacterium. Following electro-transformation, the Agrobacterium culture was incubated at 30 °C for a two-day period before being transferred to a kanamycin-containing YEB liquid medium to extend the incubation process. After an hour-long incubation at 28 °C in a shaking environment of 170 rpm, centrifugation at 4000 rpm for 4 min removed the supernatant. The organisms were then resuspended in a solution comprising 10 mM MgCl₂ and 120 μM acetylsulfate (AS), and the OD600 value of the bacterial solution was adjusted to approximately 0.6. Healthy N. tabacum plants were selected for infiltration with the bacterial solution through the lower epidermis of the leaves, using a needleless 1 mL syringe. Post injection, the N. tabacum plants were cultivated in low-light conditions for two days. Subsequently, slides were prepared from the leave's injection site to detect the GFP fluorescence signal using a laser confocal microscope (Nikon C2-ER Confocal Microscope, Tokyo, Japan).

Physiological analysis

Three A. sativa plants per treatment group were randomly selected for phenotypic observations to evaluate the effects of different treatments and time points on growth. All observations were conducted under uniform light and background conditions to ensure comparability of results. Images of the A. sativa plants were captured from both front and side angles to document morphological changes. Additionally, six A. sativa plants from each treatment group were randomly selected, their fresh leaf tissues collected and then promptly frozen in liquid nitrogen for the subsequent analysis. The GSH content was determined in shoots and roots using a GSH content assay kit (Solarbio, China). The procedure involved homogenizing 0.1 g of A. sativa tissue in 1 mL of Reagent 1, disrupting the tissue with a TissueLyser II (QIAGEN, Germany), followed by the addition of multiple reagents, with absorbance measured at 412 nm.

Prediction of upstream transcription factors and interacting proteins

Transcriptome analysis of A. sativa roots under drought and cadmium stress was performed utilizing the VennDiagram and UpSetR packages in R to create Venn diagrams and identify genes common to both differential gene sets. The motif information for transcription factor binding was then assessed from the JASPAR database (http://jaspar.genereg.net) for these shared differential genes and predicted transcriptional target genes using MEME FIMO software (version 5.5.6). Subsequently, we mapped the network of potential transcription factor interactions associated with these target genes. Additionally, we analyzed the target protein interaction networks through the STRING protein interactions database (http://string-db.org). The BLASTX was employed to compare the sequences in the target gene set against the protein sequences of reference species within the STRING database, such as A. thaliana and O. sativa, thereby constructing the interaction network based on the protein interactions of these reference species.

Statistical analysis

Figures were generated using GraphPad Prism software (version 9.0.0), with data presented as mean values and corresponding standard deviations. Statistical analyses were conducted using IBM SPSS Statistics (version 26.0). To ensure the robustness of the experimental results, three biological replicates were conducted for each treatment group in the RT-qPCR experiments, while six biological replicates were performed for each treatment group in the GSH content determination experiments. Significance of the data was assessed using one-way analysis of variance (ANOVA), with a p-value threshold of less than 0.05 established for determining statistical significance.