A mixture of Q. dentata seeds and sand with a humidity of 60% was stored in a refrigerator at 4°C for 3 weeks and then soaked in water at 40°C three times at room temperature (Zhang et al., 2023). After treatment, the Q. dentata seeds were sown in a 1:1 vermiculite-peat substrate and cultivated at 25°C/20°C with a 16 h/8 h light/dark cycle at 60% humidity. The 1/2 Hoagland nutrient solution was prepared without calcium nitrate according to the manufacturer’s protocol (Hope Bio-Technology Corporation Ltd., Qingdao, China). Seedlings of the same size and state were transplanted into 1/2 Hoagland nutrient solution and grown under the same conditions as used for sowing. Cd-grown and Pb-grown seedlings were treated with 200 mg/L CdSO₄·8/3H₂O and 1000 mg/L PbCl₂, respectively (Gao et al., 2020). Seedlings cultivated in standard 1/2 Hoagland solution were used as controls. Roots, stems, and leaves of Cd-grown, Pb-grown and control seedlings were separately harvested for RNA extraction at 0 h and 24 h after transplantation. For each treatment, three culture bottles were used, with one seedling per bottle.

Genomic data of Q. dentata were obtained from the National Genomics Data Center (https://ngdc.cncb.ac.cn/) under BioProject PRJCA013491 with accession number GWHBRAD00000000) (Wang et al., 2023). To locate QdGST genes within the genome, a BLAST search was conducted with GST genes from Arabidopsis thaliana as the query sequences. In addition, Arabidopsis GST protein sequences were retrieved from the TAIR database (https://www.arabidopsis.org/) and used to screen for Q. dentata GSTs on the basis of sequence similarity (BLASTP; e-value ≤ 10^–10). Hidden Markov model profiles for the GST_N (PF02798) and GST_C (PF00043) domains were obtained from the Pfam database and used to identify GST domains. Q. dentata protein sequences identified in homology searches were subsequently analyzed using the NCBI conserved domain database, HMMER, Pfam, and SMART (Letunic et al., 2004). The resulting QdGSTs were classified on the basis of their homology to Arabidopsis GSTs using a previously described standard method (Ghangal et al., 2020).

Sequence alignment of GST proteins from Arabidopsis and Q. dentata was performed using the MUSCLE wrapper in TBtools with default parameters to examine evolutionary relationships among the proteins (Chen et al., 2020). A maximum likelihood phylogenetic tree was built using the ‘One Step Build a ML Tree’ feature in TBtools, applying the JTT substitution model with a 95% site coverage cut-off. Node confidence was evaluated using 1000 bootstrap replicates. GST classes were visualized using with distinct colors for clarity.

The characteristics of QdGST proteins, including molecular weight, isoelectric point, and grand average of hydropathicity, were analyzed using the ExPasy tool (http://web.expasy.org/). Their subcellular localizations were predicted with Cell-PLoc 2.0 (http://www.csbio.sjtu.edu.cn/bioinf/Cell-PLoc-2/).

The chromosomal locations of 86 QdGST genes were obtained from a genome annotation file and mapped to chromosomes using TBtools (Chen et al., 2020). Genomic data for A. thaliana, C. annuum, Vitis vinifera, Triticum aestivum, and Oryza sativa were obtained from EnsemblPlants (http://plants.ensembl.org); those for Quercus mongolica were obtained from the NCBI Sequence Read Archive under accession codes PRJNA609556 and PRJNA607679. All the genomic data were analyzed for collinear relationships using the ‘One Step MCScanX’ feature in TBtools with default settings (Chen et al., 2020).

Structural information regarding the QdGST genes was obtained from a GFF file and used to visualize the conserved domains and exon–intron organization with TBtools. Conserved motif sequences and types of QdGST genes were analyzed using MEME Suite 5.5.3 (http://meme-suite.org/tools/meme) with the following parameters: motif site distribution = any number of repetitions; and maximum motif number = 10. The structures and motif distributions of QdGST genes were grouped according to a phylogenetic tree and visualized using TBtools (Chen et al., 2020).

Promoter regions (2,000 bp upstream of the translation initiation site) were extracted from all QdGST genomic sequences using TBtools (Chen et al., 2020). Potential cis-regulatory elements within these regions were identified by searching the PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) with default settings.

To investigate the tissue-specific expression of the QdGST genes, we obtained RNA expression data for Q. dentata from the National Genomics Data Center (accession code GWHBRAD00000000). For analysis of the expression of QdGST genes under Pb and Cd stresses, samples were collected using the method described above. Total RNA extraction and mRNA library construction and sequencing were performed according to the methods described by Chen et al. (2022). All RNA-seq data have been submitted to the Genome Sequence Archive (accession CRA013085) of the National Genomics Data Center (CNCB-NGDC Members and Partners, 2022). Expression levels were normalized to FPKM (fragments per kilobase million), and a heatmap was generated by applying the logarithmic transformation log₁₀(FPKM+1) using TBtools.

Total RNA was extracted from the samples using a SteadyPure Plant RNA Extraction Kit (Accurate Biotechnology, Hunan, China) according to the protocol from manufacturer. Complementary DNA (cDNA) was synthesized using Evo M-MLV RT Premix for qPCR (Accurate Biotechnology). The primers for the qRT-PCR experiments were designed using Primer3plus (https://www.primer3plus.com/) and are listed in Supplementary Table S3 . The CaCs gene was selected as the internal reference gene for roots, whereas the EF1-α gene was used for stems and leaves ( Supplementary Table S3 ). Subsequently, the qPCRs were performed on a Bio-Rad CFX96 Touch Real-Time PCR Detection System (USA) using an SYBR Green Premix ProTaq HS qPCR Kit (Accurate Biotechnology) under the following conditions: 95°C for 30 s, followed by 40 cycles of 95°C for 15 s and 60°C for 30 s. The 2^−ΔΔCt method with three technical replicates per sample was used to calculate relative gene expression (Livak and Schmitten, 2001).

Gene-specific primers ( Supplementary Table S3 ) for QdGSTU20 and QdGSTU36 were used to amplify their coding DNA sequence regions. To enable expression in yeast (Saccharomyces cerevisiae), cDNAs for QdGSTU20 and QdGSTU36, obtained through PCR amplification, were cloned into the KpnI-XbaI restriction sites of the pYES2.0 plasmid. A transformation kit (Huayueyang Biotechnology, China) was used to introduce the empty plasmid and recombinant vector into S. cerevisiae INVSc1. Jiang’s method with minor modifications was used to evaluate the metal resistance of the transgenic yeast (Jiang et al., 2024). Yeast harboring either the recombinant plasmid or the empty plasmid was cultured in liquid medium (SC-U/Glu) and incubated at 28°C and 200 rpm until the optical density at 600 nm reached 1.0. Subsequently, serial ten-fold dilutions were performed using sterile water, and 2-μL aliquots from each dilution were dropped onto solid SC-U/Gal induction medium supplemented with various heavy metals (10 μM CdSO₄·8/3H₂O, 3 mM PbCl₂, or 15 mM MnSO₄·H₂O) or control medium (no additional metal ions). Three independent biological replicates were included in each treatment. Growth phenotypes were recorded after 3 days of incubation at 28°C.

GraphPad Prism 9.4.1 was selected for statistical analysis and graph generation. In figures showing the results of qPCR analyses, error bars represent the standard deviation from three independent biological replicates. Statistical analysis was performed using one-way analysis of variance and post hoc least significant difference tests to identify differences among group means, at a significance level of 0.05.

Eighty-six GSTs in Q. dentata were identified and systematically classified on the basis of their chromosomal locations and the homology of their encoded proteins with that of Arabidopsis GSTs ( Supplementary Table S1 ) (Zhao et al., 2021; Mo et al., 2023). GST proteins from Q. dentata (86) and Arabidopsis thaliana (56) were used to construct a phylogenetic tree for investigation of the phylogenetic relationships among the QdGSTs. QdGSTs belonging to six classes (tau, lambda, phi, TCHQD, theta, and DHAR) were identified ( Figure 1 ), with the majority belonging to the tau class (56 members, 65.1%), followed by the lambda (13 members, 15.1%), phi (eight members, 9.3%), and theta (five members, 5.8%) classes. Only two members each belonged to the TCHQD and DHAR classes, and no QdGST belonged to the zeta class ( Figure 1 ).

The coding DNA sequences of QdGSTs varied in length from 426 bp (QdGSTU39) to 1467 bp (QdGSTT1), encoding proteins that ranged from 141 to 488 amino acids. The molecular weights (MWs) of the QdGSTs ranged from 15.9 kDa (QdGSTU39) to 55.1 kDa (QdGSTT1), and their isoelectric point values ranged from 4.48 (QdGSTL1) to 9.38 (QdTCHQD1 and QdTCHQD2). The grand average of hydropathicity values of all QdGSTs, except for QdGSTU3 (0) and QdGSTU54 (0.133), were negative and ranged from −0.02 (QdGSTU24) to −0.57 (QdMTP8.3). Most of the QdMTPs had weak hydrophilicity, and the majority were predicted to have cytoplasmic localization, with limited numbers in the chloroplast and the nucleus ( Supplementary Table S1 ).

The 86 QdGST genes were unevenly dispersed across 11 of the 12 chromosomes. Twenty-one (24.4%) and 20 (23.3%) genes, belonging to the tau and lambda members, were located on chromosomes 4 and 7, respectively; chromosome 9 contained 13 QdGST genes; and chromosomes 6, 2, 3, 12, 10, and 8 had 8, 7, 6, 4, 3, and 2 QdGST genes, respectively. Chromosomes 5 and 11 contained only one QdGST gene each ( Figure 2 ). No significant association was observed between chromosomal length and number of QdGST genes ( Figure 2 ). Segmental and tandem duplication have important roles in driving the expansion of gene families (Wang et al., 2019). A total of 31 gene pairs distributed on seven chromosomes were identified as tandem duplication types; 22 of these (71.0%) involved genes from the tau class, and seven (22.6%) involved genes from the lambda class. There was one tandem duplication (0.03%) corresponding to each of the theta and phi classes ( Figure 2 ). In addition, seven QdGST gene pairs (QdGSTU6/QdGSTU54, QdGSTU31/QdGSTU48, QdGSTU34/QdGSTU49, QdGSTF3/QdGSTF5, QdGSTL3/QdGSTL13, QdGSTL7/QdDHAR2, and QdGSTL7/QdGSTL13), located on chromosomes 3, 6, 7, 8, 9, 10, and 12, respectively, were identified as segmental duplications ( Figure 3 ). Three of these gene pairs were from the tau class, two from the lambda class, and one from the phi class. Notably, one segmental duplication pair (QdGSTL7/QdDHAR2) contained genes from both the lambda and DHAR subfamilies. QdGSTU31, QdGSTU34, and QdGSTL7 underwent both segmental and tandem duplication. Overall, these results suggest that the duplication events identified here contributed primarily to the expansion of the tau and lambda classes of GST genes in Q. dentata.

To explore the evolutionary relationships of GST genes in Q. dentata and other species, we constructed a synteny map based on six species (A. thaliana, O. sativa, C. annuum, T. aestivum, V. vinifera, and Q. mongolica) and previously published genome-wide GST data (Wang et al., 2019; Islam et al., 2019; Li et al., 2022). We identified 135 orthologous GST gene pairs between Q. dentata and the six angiosperm species ( Figure 4 ). Of these, 42 collinear blocks of QdGST genes were found between Q. dentata and Q. mongolica, followed by 38 between Q. dentata and V. vinifera, 23 between Q. dentata and C. annuum, and 20 between Q. dentata and A. thaliana. In addition, eight and four orthologous GST gene pairs were detected for T. aestivum and O. sativa, respectively ( Figure 4 ). No QdGST gene syntenic regions were found on chromosome 5 or 11 of the Q. dentata genome. Three QdGST (QdGSTL13, QdGSTL7, and QdGSTU46) gene pairs were detected between Q. dentata and two monocot plants (T. aestivum and O. sativa); this indicates that the majority of QdGST genes were formed after the divergence of their common ancestor. Notably, 36 QdGST genes lacked detectable orthologs across other species.

We analyzed the composition of QdGST proteins and found ten conserved motifs among the 86 QdGSTs. The presence of these motifs varied among QdGST from different classes, with motif 9 found only in the lambda class, whereas motifs 5 and 10 were specific to the tau class. Motifs 4, 1, and 3 were consistently present in all proteins of the DHAR and TCHQD classes; all members of the theta class had motifs 1 and 3; and all members of the phi class had motifs 4, 1, 2, and 3. All members of the lambda class except for QdGSTL1 harbored motifs 9, 4, 1, 2, and 7; and most members of the tau class (50/56, 89.3%) contained motifs 4, 6, 1, 2, 8, and 3 ( Figure 5A ). To investigate the evolution of the QdGST genes, we analyzed their structures using the Q. dentata genome annotation file. We found that the number of introns in QdGST genes varied from 0 to 12, with genes within the same class exhibiting similar structural patterns. Most QdGST genes (49/56) in the tau subfamily had one intron, and five genes had two introns, and the QdGSTU2 gene had 11 introns. Most members (11/13) of the lambda class had 9–10 introns. Seven of eight QdGST genes in the phi class had two introns; the exception was QdGSTF8, which had three. TCHQD- and DHAR-class genes contained three and five introns, respectively, and intron numbers in the theta class ranged from four (QdGSTT1) to 12 (QdGSTT5) ( Figure 5B ).

We screened the 2,000-bp promoter regions upstream of the translation initiation sites of the 86 QdGST genes using PlantCARE to detect cis-acting elements. The results revealed 1,924 cis-regulatory elements that could be classified into 29 categories ( Figure 5C ) on the basis of their functions. These functions included response to light, plant hormones, stresses, and plant development. More than 50% of the identified elements were related to light response; these included ACE, G-box, 4cl-CMA2b, MRE, Box 4, ATCT-motif, ATC-motif, chs-CMA2a, AT1-motif, AE-box, and ACA-motif, which were ubiquitously present in the promoter regions of all QdGSTs ( Figure 5C , Supplementary Table S2 ).

A total of 527 cis-acting elements were found to be related to signaling by plant hormones, including gibberellin (TATC-box, P-box, GARE-motif), salicylic acid (TCA-element), abscisic acid (ABRE), auxin (TGA-element, AuxRR-core, AuxRE, and TGA-box), and methyl jasmonate (TGACG-motif and CGTCA-motif), as well as zein metabolism regulation (O2-site). A total of 326 cis-regulatory elements were implicated in stress responses such as anaerobic induction (ARE, GC-motif), drought response (MBS), low-temperature response (LTR), defense and stress response (TC-rich repeats), and heat response (HSE). Twenty-four and 83 cis-acting elements were related to circadian control and plant development (CAT-box, GCN4_motif, RY-element, and HD-Zip 1), respectively ( Figure 5C , Supplementary Table S2 ).

The expression patterns of QdGST genes were characterized using transcriptome sequencing data from diverse tissues and developmental stages. We detected expression (log₁₀(FPKM+1) > 0) of all QdGST genes except QdGSTU41 in at least one tissue at each development stage ( Figure 6A ). Nine QdGST genes (QdDHAR2, QdGSTL3, QdGSTL7, QdGSTL12, QdGSTL6, QdGSTU5, QdGSTU11, QdGSTT2, and QdGSTT3) showed high expression in all samples (log₁₀(FPKM+1) > 1). The eighty-six QdGST genes could be divided into three classes according to their tissue-specific expression profiles: (a) QdGST genes with extremely low expression in almost all tissues; (b) QdGST genes exhibiting low-to-medium expression across different tissues; and (c) QdGST genes with high expression in some tissues. These findings were consistent with previous reports (Islam et al., 2019; Hasan et al., 2021).

To investigate the responses of QdGST genes to heavy metal stresses, we analyzed their expression profiles under Pb and Cd treatment using data from transcriptome sequencing. Nine QdGST genes of the tau class (QdGSTU19, QdGSTU20, QdGSTU25, QdGSTU30, QdGSTU31, QdGSTU32, QdGSTU34, QdGSTU42, and QdGSTU43) and one member of the lambda class (QdGSTL12) were upregulated in response to both Pb and Cd treatment. Two members of the tau class (QdGSTU18 and QdGSTU40), one of the theta class (QdGSTT4), and one of the phi class (QdGSTF8) were upregulated only in response to Pb treatments, whereas 13 tau-class (QdGSTU2, QdGSTU17, QdGSTU22, QdGSTU29, QdGSTU35, QdGSTU36, QdGSTU37, QdGSTU41, QdGSTU44, QdGSTU47, QdGSTU50, QdGSTU55, and QdGSTU56) and four lambda-class (QdGSTL3, QdGSTL6, QdGSTL7, and QdGSTL11) genes were upregulated only by Cd treatment. QdGSTU19, QdGSTU20, QdGSTU36, QdGSTU39, QdGSTU44, QdGSTU45, and QdGSTL11 showed particularly marked upregulation (>100-fold increase in expression) under Pb and/or Cd treatment ( Figure 6B ).

The results of the RNA-seq analysis for eight QdGST genes with high expression levels under Pb and Cd treatments were validated using qRT-PCR. Under normal growth conditions, the expression patterns of these genes showed tissue specificity, with QdGSTL6, QdGSTL12, QdGSTU17, and QdGSTU36 having higher levels in stems and leaves than in roots ( Figure 7 ), whereas QdGSTU20 and QdGSTU44 had higher expression levels in roots and stems compared to leaves, and QdGSTU19 and QdGSTU34 were highly expressed in stems and roots, respectively. All the tested QdGST genes were upregulated under both Pb and Cd treatments in stems, and all except QdGSTL6 were upregulated under Cd treatment in roots. Expression of QdGSTL12, QdGSTU19, QdGSTU20, QdGSTU34, QdGSTU36 and QdGSTU44 was induced by Pb treatment. In leaves, the expression of QdGSTU20, QdGSTU34, QdGSTU36, and QdGSTU44 was induced by Cd treatment, whereas that of QdGSTL6 was upregulated by Pb treatment. However, four QdGST genes (QdGSTL6, QdGSTL12, QdGSTU17, and QdGSTU19) were inhibited by Cd and three (QdGSTL12, QdGSTU20, and QdGSTU34) by Pb in leaves.

To further characterize the involvement of QdGST genes in heavy metal resistance, QdGSTU20 and QdGSTU36, which were highly expressed under both Pb and Cd treatments, were introduced into wild-type yeast for metal sensitivity analysis. Under normal conditions, yeast strains with empty vectors and those expressing QdGSTU20 or QdGSTU36 showed no significant difference in growth. However, when expressed in yeast, QdGSTU36 conferred Pb and Cd tolerance. Expression of QdGSTU20 did not alter the sensitivity of yeast to Pb or Cd; however, it did confer Mn tolerance ( Figure 8 ).