Sweet potato (I. batatas L.) cv. ‘Haida HD7791’ was used in this study, following the hydroponic protocol described by (Kumar et al., 2024). Briefly, stem cuttings were surface-sterilized with Carbendazim (1 g L^-1), rooted in distilled water, and then transferred to half-strength Hoagland nutrient solution. The nutrient solution pH was monitored daily with a calibrated pH meter and adjusted to 6.0 ± 0.5 with dilute NaOH or HCl in a controlled growth room (25 to 27°C, 16 h light/8 h dark). After 6 days, seedlings were subjected to pre-treatments with MT (1 µM) and GSH (2 mM) for 7 days, followed by exposure to 40 μM Cr(VI) applied as K₂Cr₂O₇ in the same nutrient solution. The RNA-seq experiment used four treatments: CK, Cr, MT + Cr (MC), and GSH + Cr (GC), each with three biological replicates of six plants per replicate. Growth, physiological and biochemical traits, including plant height, leaf number and area, shoot and root fresh/dry weight, relative water content (RWC), root architectural traits, gas-exchange parameters, chlorophyll fluorescence and pigments, ROS and MDA, electrolyte leakage, osmoprotectants (soluble sugars and proline), antioxidant metabolites (GSH, GSSG, AsA, phenolics, and flavonoids), antioxidant enzymes, stomatal traits, and Cr uptake/translocation, were measured as previously described by (Kumar et al., 2024), which were subsequently used as external traits in the WGCNA analysis.

For RNA-seq, leaf tissue was collected from plants grown under the four treatments (CK, Cr, MC, and GC) at the same developmental stage used for physiological measurements. For each treatment, samples from multiple plants within a replicate pot were pooled to form one biological replicate, yielding 12 libraries in total. Total RNA was extracted using a plant RNA isolation kit (RNA Easy Fast Plant Tissue, Tiangen Biotech, Beijing, China) according to the manufacturer’s instructions. RNA integrity and purity were evaluated by 1% agarose gel electrophoresis, spectrophotometry (A260/280 and A260/230 ratios), and an Agilent 2100 Bioanalyzer. For each qualified sample, mRNA was purified from 1-2 µg of total RNA using oligo(dT) magnetic beads, fragmented, and reverse transcribed to first-strand cDNA using random hexamer primers, followed by second-strand synthesis. After end repair, A-tailing, adapter ligation, and PCR enrichment, size-selected cDNA fragments were used to construct strand-specific libraries following the standard Illumina protocol. The libraries were quantified and sequenced on an Illumina HiSeq X Ten platform (Illumina, San Diego, CA, USA) to generate paired-end reads.

FastQC v0.11.9 (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) was used to assess raw read quality, including per-base quality scores, GC content, sequence length distribution, and adapter contamination. Raw reads were then filtered with Trimmomatic v0.39 (Bolger et al., 2014) to remove Illumina adapter sequences, trim low-quality bases, discard reads containing > 10% unknown bases (N), and eliminate reads shorter than 36 bp. Reads were reevaluated via FastQC to confirm the removal of low-quality and contaminated sequences, and only high-quality clean reads were retained for downstream analyses. Clean reads from each library were aligned to the sweet potato reference genome Xushu 18 genome with the annotation v1.0.a2 using HISAT2 (v2.2.1) with default parameters (Kim et al., 2019), and the resulting SAM files were converted, sorted, and indexed as BAM files with SAMtools v1.10 (Danecek et al., 2021). Gene-level read counts were obtained with featureCounts v2.0.0 (Liao et al., 2014) using the corresponding genome annotation, while transcript-level abundance was estimated as transcripts per million (TPM) with StringTie v2.2.1 (Pertea et al., 2015) to provide normalized expression values for subsequent differential expression and co-expression analyses. RNA-seq data has been submitted to the National Genomics Data Center (NGDC) database (accession number: PRJCA053307; https://ngdc.cncb.ac.cn).

Genes with very low mean normalized expression were removed, and 68,922 genes were retained for downstream analysis. Principal component analysis (PCA) and hierarchical clustering were performed in R 4.4.1 (http://www.r-project.org/) to assess sample variability and detect potential outliers. Differential expression analysis for all pairwise contrasts: T1 (Cr vs Ck), T2 (GC vs Ck), T3 (MC vs Ck), T4 (GC vs Cr), T5 (MC vs Cr), and T6 (GC vs MC), was carried out with DESeq2 (Love et al., 2014) based on raw gene-level counts. P-values were adjusted using the Benjamini–Hochberg method, and genes with |log2FC| ≥ 1 and FDR ≤ 0.05 were defined as differentially expressed. Heatmaps and volcano plots were generated using the R packages pheatmap and ggplot2, respectively. For expression clustering, all DEGs were clustered using visCluster in the ClusterGVis package v0.1.0 (https://github.com/junjunlab/ClusterGVis), with normalized TPM values as input, to identify major expression patterns across CK, Cr, GC, and MC treatments. Transcription factors among DEGs and hub genes were annotated using PlantTFDB (Jin et al., 2017).

All expressed genes and DEGs were functionally annotated against public databases, including NCBI non-redundant (NR), Swiss-Prot, Gene Ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG), using BLAST-based workflows and eggNOG-mapper. Functional enrichment analysis of DEGs, co-expression clusters, and trait-related WGCNA modules was performed in R 4.4.1 with clusterProfiler (Wu et al., 2021), considering GO biological process (BP), molecular function (MF), cellular component (CC) terms, and KEGG pathways (Kanehisa and Goto, 2000). The full set of expressed genes served as the background, and the Benjamini–Hochberg correction was applied to control the FDR; GO terms/pathways with adjusted p (FDR) ≤ 0.05 and at least 2 genes per term were considered significantly enriched. The top enriched categories were visualized using ggplot2 and enrichplot.

Weighted gene co-expression network analysis (WGCNA) was performed in R 4.4.1 using the WGCNA package v1.73 (Langfelder and Horvath, 2008) to identify modules associated with Cr tolerance and glutathione-mediated physio-chemical improvements. Genes with extremely low expression or variance were removed, and the remaining expression matrix (selected top 25%) across all samples was used as input. A signed network was constructed by calculating pairwise Pearson correlation coefficients between genes and transforming them into an adjacency matrix using a soft-thresholding power selected with the pickSoftThreshold function to approximate scale-free topology (R² close to 0.8). The adjacency matrix was then converted into a topological overlap matrix (TOM), and the corresponding dissimilarity (1–TOM) was used for average-linkage hierarchical clustering. Gene modules were defined using the dynamic tree cut algorithm with a minimum module size of 30 genes, and closely related modules (eigengene dissimilarity < 0.25) were merged. Module eigengenes (MEs) were correlated with the panel of growth, physiological, and biochemical traits. Modules significantly associated with improved photosynthetic and antioxidant traits and/or reduced oxidative damage and Cr load were considered GC-related modules. Within these modules, key genes were defined based on high module membership (MM) and gene significance (GS) values (e.g., |MM| ≥ 0.8 and |GS| ≥ 0.3). The significant modules were exported to Cytoscape (Otasek et al., 2019), and the CytoHubba plugin (Chin et al., 2014) was used to calculate node degree and identify the top 30 hub genes, which were further analyzed as candidate regulators underlying glutathione-driven physio-chemical improvements under Cr stress.

To validate the RNA-seq results, eight potential candidate genes were selected for RT-qPCR analysis. Leaf samples from the four treatments mentioned above were collected in parallel with those used for RNA-seq, with three biological replicates per treatment. Total RNA was extracted as described above, and 1 µg of DNase-treated RNA was reverse transcribed using a commercially available cDNA synthesis kit (Fastking gDNA Dispelling RT SuperMix, Tiangen Biotech, Beijing, China) according to the manufacturer’s instructions. Gene-specific primers were designed in Primer3 software, and a housekeeping gene (Actin) was used as an internal reference. qRT-PCR was performed using a MA-6000 PCR system (Suzhou Yarui Biotechnology, China). Each sample was run in triplicate. Relative expression levels were calculated using the 2^-ΔΔCt method. The primers of selected genes are listed in Supplementary Table 1.

We sequenced 12 RNA-Seq libraries to decode the molecular crosstalk between GSH and MT to mitigate the effect of Cr stress. Across 12 libraries, sequencing quality was high, with 99.60% to 99.70% of raw reads retained as clean reads (Supplementary Table 2). The overall genome alignment was consistent across treatments, with total mapped reads ranging from ~84% to ~86%, of which about 64% to 66% were uniquely mapped and ~20% mapped to multiple loci. Unmapped reads remained low (14% to 16%), indicating reliable RNA-seq data suitable for downstream differential expression analysis. Given the high genome complexity of sweet potato (a highly heterozygous hexaploid, 2n = 6x = 90), mapping rates of 84–86% are expected, as polyploidy and repetitive/homologous sequences can reduce unique alignment. PCA accounted for 88.4% and 10.3% of the total variance along PC1 and PC2, respectively (Figure 1a). The CK clustered together, confirming replicate consistency. Cr stress induced a substantial transcriptional shift relative to CK, while MC clustered near Cr, indicating limited transcriptomic deviation from Cr stress. In contrast, GC samples were clearly separated from Cr along PC1, demonstrating a distinct expression pattern. These results suggest that GSH exerts a more pronounced effect on transcriptomic reprogramming during Cr detoxification than melatonin. Hierarchical clustering of transcriptomes revealed a clear treatment-level structure (Figure 1b), with biological replicates clustering tightly within each condition, confirming reproducibility. GC formed a distinct branch, indicating a transcriptional program that diverges from both CK and Cr stress. In contrast, MC grouped with Cr on the same major branch, consistent with a MT response that remains close to the Cr-stress state. Pairwise Pearson correlations supported this hypothesis (Figure 1c). MC vs Cr showed a robust correlation (r = 0.97), indicating minimal deviation from chromium-induced expression. GC vs Cr correlations were lower (r = 0.92–0.94), indicating a larger GC-driven reprogramming of the Cr response. GC vs MC correlations were lower still (r = 0.87–0.89), highlighting distinct effects of GSH and MT under identical Cr exposure. Compared with CK, GC showed a lower positive correlation than either Cr or MC, consistent with GC entering a unique detoxification state rather than reverting to a control-like profile. The boxplots of all groups showed variation among the samples and treatment groups (Figure 1d). Mean expression comparisons were significant for all pairwise contrasts (p ≤ 0.05) except MC vs Cr, which was non-significant, again indicating that MT-treated leaves were close to Cr stress, whereas GSH drives a more distinct adjustment. These results demonstrated that GC has a stronger, more distinct expression profile than Cr and MC, partially supporting our hypothesis.

Using a threshold of log₂FC ≥ 1 and FDR ≤ 0.05, Cr stress versus CK (T1) yielded 5,601 DEGs, of which 56.60% were upregulated, and 43.40% were downregulated (Figures 2a, b; Supplementary Table 3). In plants applied with GSH under Cr (T2: GC vs CK), 4,677 DEGs were detected, with a nearly balanced distribution of 49.10% upregulated and 50.90% downregulated genes. MT under Cr (T3: MC vs CK) produced 5,401 DEGs, with 55.60% upregulated and 44.4% downregulated, indicating a transcriptional response similar in magnitude to Cr alone (Figures 2a, b). Direct comparisons with Cr alone revealed that GSH (T4: GC vs Cr) reprogrammed 824 genes, with a strong bias toward repression (17.50% up- vs 82.50% downregulated), whereas MT (T5: MC vs Cr) affected only 319 genes (25.10% up- and 74.90% downregulated). The GC vs MC contrast (T6) identified 533 DEGs, of which 64.50% were higher in GC, and 35.50% were lower. These DEG patterns are consistent with our hypothesis that although both MT and GSH act in a Cr-stress context, GSH exerts a quantitatively stronger and more distinct modulation of the Cr-responsive transcriptome than MT (Figures 2a–c). The top upregulated and downregulated genes showed different expression profiles under different groups (Figures 2a–c).

GO enrichment analysis of DEGs (FDR ≤ 0.05) for the three treatments (T1 to T3) versus CK showed a common pattern of activation of detoxification/transport processes and repression of photosynthetic and growth-related functions (Figure 3; Supplementary Figure S1; Supplementary Table 4). In T1 (Cr vs CK), upregulated genes were mainly involved in xenobiotic transmembrane transport, toxin and secondary-metabolite metabolism, responses to wounding, fungus and abscisic acid, monocarboxylic acid transport, and cellular nutrient responses (Figure 3a). These genes were enriched in peroxisomes, chloroplasts, and other membranes, and ABC transporter complexes, with molecular functions dominated by various transmembrane transporters and oxidoreductase activities. Downregulated genes in T1 were strongly enriched in steroid/terpenoid biosynthesis and metabolism, regulation of hormone levels, cell-wall macromolecule metabolism/organization, and photosynthesis and light reactions, together with chloroplast and thylakoid membrane terms and peroxidase/glucosyltransferase activities (Figure 3b). Likewise, in T2 (GC vs CK), the upregulated DEGs showed a similar but more detoxification-focused profile, with significant enrichment of response to ROS, cellular aldehyde metabolic process, detoxification, terpenoid metabolism, ABA response, monocarboxylic acid transport/catabolism, and polyol metabolism (Figure 3c). In contrast, downregulated GC-responsive genes were enriched in steroid/terpenoid and glucan/oligosaccharide biosynthesis, cell-wall organization, photosynthesis, and photosynthetic membrane components, indicating that photosynthetic machinery remains lower than in CK despite GSH treatment (Figure 3d). In T3 (MC vs CK), upregulated genes were enriched for xenobiotic metabolism and transport, toxin metabolism, and transporter/oxidoreductase activities associated with peroxisomes, and ABC transporter complexes (Supplementary Figure 2). In contrast, downregulated genes showed the same core categories as T1 and T2; steroid/terpenoid biosynthesis, photosynthesis and photosynthetic membranes, and glucan/cell-wall metabolism (Supplementary Figure 2). These results indicate that the shared GO signatures confirm that Cr stress suppresses chloroplast and primary metabolic functions. However, GSH particularly strengthens ROS and aldehyde detoxification pathways compared with Cr alone and MT (Figure 3; Supplementary Figure S2).

In the direct comparisons under Cr stress, GO enrichment highlighted that GSH and MT differentially remodel the stress transcriptome (Supplementary Figures 3, 4). Relative to Cr alone (T4), GC upregulated genes associated with peptidase regulation and inhibitor activity, as well as acetyl-/acyltransferases and glyoxalase enzymes, indicating enhanced protein quality control and methylglyoxal/ROS detoxification. Conversely, GC downregulated genes related to cell-wall formation, apoplast/cell-surface signaling, and solute symport, consistent with a reduced need for cell-wall remodeling and stress signaling once Cr toxicity is alleviated. Compared with Cr (T5), MC mainly induced genes involved in membrane and vesicle transport and metal/amine ion transport, while repressing amino-acid, nitrogen, and sulfur compound biosynthesis. This suggests that MT adjusts intracellular trafficking and secondary metabolism but invests less in sulfur-based detox pathways than GC. The GC vs MC contrast (T6) further showed that GC preferentially upregulated ion and anion homeostasis, including nitrite efflux and inorganic anion transport. In contrast, genes linked to carbohydrate/biogenic amine metabolism and proteasome/peptidase regulation were relatively higher in MC. These GO patterns support the view that GSH more strongly reinforces redox and detoxification networks and ion homeostasis than MT under the same Cr stress (Supplementary Figures 2, 3).

KEGG enrichment of DEGs (FDR ≤ 0.05) for T1-T3 versus CK showed that Cr stress mainly activated stress-responsive lipid and secondary metabolism while suppressing photosynthetic carbon metabolism (Figure 4; Supplementary Figure S3; Supplementary Table 5). In T1 (Cr vs CK), upregulated genes were enriched in sesquiterpenoid and triterpenoid biosynthesis, α-linolenic acid metabolism, zeatin and other secondary-metabolite biosynthesis, MAPK signaling, glutathione metabolism, and peroxisome pathways, consistent with enhanced hormone/oxylipin signaling and ROS detoxification. Conversely, Cr strongly downregulated photosynthesis and antenna proteins, starch and sucrose metabolism, carbon fixation, porphyrin and carotenoid biosynthesis, and several sugar and steroid biosynthetic pathways, reflecting inhibition of light harvesting and primary carbon metabolism. In T2 (GC vs CK) and T3 (MC vs CK), upregulated pathways largely overlapped with the Cr response, including terpenoid and flavonoid/phenylpropanoid biosynthesis, α-linolenic acid metabolism, fatty-acid degradation, and MAPK signaling. However, GC uniquely enriched cysteine and methionine metabolism, glycolysis/gluconeogenesis, pyruvate and nitrogen metabolism, and protein processing in the ER, together with glutathione metabolism, indicating that exogenous GSH more strongly reinforces sulfur assimilation, central carbon metabolism, and the GSH pool required for detoxification.

In the direct comparisons under Cr stress, KEGG enrichment highlighted apparent differences between GSH and MT responses (Supplementary Figures 3, 4). Relative to Cr alone (T4), GSH-upregulated genes were mainly enriched in glutathione metabolism, cysteine and methionine metabolism, and related carbon/nitrogen metabolism pathways, consistent with reinforced sulfur-based antioxidant capacity and detoxification. Pathways reduced by GC included several cell wall- and stress-signaling-related routes, in line with a reduced stress load once Cr toxicity is alleviated. In contrast, MT vs Cr (T5) mainly enriched pathways associated with plant hormone signal transduction, MAPK signaling, and membrane/vesicle transport, suggesting that MT primarily fine-tunes signaling and trafficking rather than strongly boosting core detox metabolism. The GC vs MC contrast (T6) further showed that genes higher in GC were enriched in glutathione and sulfur metabolism, nitrogen metabolism, and other redox-related pathways. In contrast, genes higher in MC were more associated with secondary metabolism and protein processing.

To obtain a mechanistic view of the transcriptomic response, we clustered all DEGs across the six comparisons, yielding 7,734 genes that were grouped into four major expression clusters (C1–C4) (Figure 5; Supplementary Table 6). Cluster sizes ranged from 1,234–3,285 genes, with C1 containing 1,521 genes (19.7%), C3 containing 1,694 genes (21.9%), C4 containing 1,234 genes (16.0%), and C2 containing the largest group with 3,285 genes (42.5%). The z-score heatmap and average expression profiles clearly separated Cr stress (CR) from the control (CK) and revealed distinct GC- and MC-modulated patterns under the same Cr background. C1 comprised genes strongly induced by Cr relative to CK, whose expression was partially reversed by GC and, to a slight extent, by MC. GC brought these Cr-upregulated genes back to basal levels, whereas MC maintained a moderate induction, indicating that GC suppresses a larger subset of Cr-activated transcripts more efficiently. C3 also contained Cr-induced genes, but MC further enhanced their expression beyond Cr, while GC slightly reduced their levels. This pattern defines an MC-preferential module, suggesting that MT tends to amplify specific Cr-responsive pathways rather than reduce them. C4 represented another Cr-induced group in which GC further boosted expression above Cr, with MC showing intermediate induction. These genes likely correspond to GC-enhanced detoxification and defense programs that are selectively activated in the presence of GSH. In contrast, C2 genes were highly expressed in CK but strongly repressed by Cr, and remained low under both GC and MC treatments, indicating a large cohort of Cr-suppressed, putatively growth-associated genes that neither treatment fully restores. The clustering analysis shows that, although both GC and MC act on Cr-responsive genes, GSH predominantly counteracts harmful Cr-induced expression (C1) and reinforces a specific set of protective genes (C4). In contrast, MT more often maintains or even strengthens Cr-like expression states (C3). These coordinated expression modules provide a mechanistic framework supporting GSH as the more effective regulator of the Cr-stress.

Consistent with these expression patterns, GO and KEGG enrichment confirmed that each cluster represents a functionally coherent module (Figure 5). C1 was enriched for stress- and defense-related GO terms such as response to wounding, chitin, bacteria, and metal ions, together with KEGG pathways including α-linolenic acid metabolism, sesquiterpenoid/triterpenoid, phenylpropanoid, and flavonoid biosynthesis, and MAPK signaling pathway, indicating a broad Cr-induced defense program that is reduced mainly by GC. C2 was strongly enriched for photosystem and chloroplast thylakoid components and for photosynthesis-related processes, with KEGG terms including photosynthesis, antenna proteins, starch and sucrose metabolism, carbon fixation, and porphyrin metabolism, consistent with a core photosynthetic module that is repressed by Cr and not fully restored by either treatment. C3 showed GO enrichment for xenobiotic transport and toxin/xenobiotic metabolism, and KEGG enrichment for fatty-acid degradation, α-linolenic acid and sphingolipid metabolism, nitrogen metabolism, and peroxisome, defining an MC-preferential detoxification and lipid–nitrogen remodeling module. Finally, C4 was enriched for responses to ROS, H₂O₂, heat, and chloroplast organization, with KEGG pathways such as starch and sucrose metabolism, glycerophospholipid and pyruvate metabolism, secondary-metabolite biosynthesis, and protein processing in the ER, highlighting a GC-enhanced module that combines redox detoxification with metabolic and protein-quality control adjustments under Cr stress.

To further dissect the coordinated transcriptional responses underlying Cr detoxification by GSH and MT, we constructed a co-expression network using the top 25% most variable genes. A soft-thresholding power of β = 12 produced a scale-free topology with R² = 0.90, indicating a biologically meaningful network structure (Figure 6a). WGCNA grouped the genes into seven co-expression modules with dynamic tree-cut (blue, yellow, black, pink, turquoise, brown, green, and grey), and their eigengenes were correlated with growth, photosynthetic, oxidative-stress, and Cr-accumulation traits (Figure 6b). Among these, the blue module showed the most distinct pattern and its eigengene was strongly and positively correlated with biomass and water-status traits (SFW, SDW, RFW, RDW, and RWC; r = 0.72-0.99, p ≤ 0.01) and with photosynthetic performance (Pn, Gs, Ci, Trmmol, Y(II), qP, ETR, and total chlorophyll; r = 0.76-0.99, generally p < 10^-3) (Figures 6c, d). In contrast, the blue eigengene was negatively correlated with membrane injury and oxidative markers (EL, O₂^-, MDA, and H₂O₂; r = -0.84 to -0.90, p ≤ 10^-3) and showed robust negative correlations with Cr accumulation and translocation parameters (LU, SU, RU, and Cr_trans; r = −0.98 to −1.00, p ≤ 10^-7). Modules such as yellow and turquoise displayed an opposite trend, as they were negatively associated with growth and photosynthetic traits and positively associated with oxidative damage and Cr accumulation, suggesting that these modules represent Cr-sensitive expression programs. These results highlighted that the blue module is a core Cr-tolerance module and upregulation is associated with improved growth, photosynthetic efficiency, and reduced Cr uptake/translocation. Therefore, these physiological improvements are most pronounced under GSH treatment; the blue module likely captures a significant portion of the GC-driven protective network against Cr stress.

Following WGCNA analysis, we selected the top 30 hub genes by degree as key candidates underlying the physio-chemical improvements observed in GSH-treated plants (Figure 7; Supplementary Table 7). Several hub genes, such as HEMA1, CHLP, CAB4, CAB13, PSAN, and PSAL, encode enzymes and structural proteins involved in chlorophyll/tetrapyrrole biosynthesis and photosystem I light harvesting, directly linking the GC module to the recovery of chlorophyll content, Fv/Fm, and electron transport efficiency under Cr stress. A second group of hubs is closely tied to redox and detoxification processes that cooperate with glutathione metabolism, including a peroxidase annotated with “response to oxidative stress,” “hydrogen peroxide catabolic process,” and phenylpropanoid biosynthesis, and CHLP/HEMA1 with GO terms related to cell redox homeostasis and porphyrin metabolism. These genes likely work in concert with the AsA-GSH cycle to remove ROS and toxic compounds in GC-treated plants. Transport- and osmotic-regulation hub genes, such as aquaporin TIP2-1 (two copies) and the phosphate transporter PHT1-4, together with TPPA (trehalose-phosphate phosphatase), are associated with improved water status, nutrient uptake, and osmoprotectant (trehalose) metabolism, consistent with higher RWC, better root architecture, and reduced EL. Growth, signaling, and structural adjustment are represented by hub genes, such as GA20OX1/GA20OX2 (gibberellin biosynthesis), the HD-ZIP transcription factor HAT4 (light and hormone-responsive developmental regulator), a receptor-like serine/threonine kinase, and multiple cell-wall and lipid–related enzymes (XTH1, pectate lyases, alpha-mannosidase, ASAT1, and CAS1). These genes integrate hormone signaling, cell-wall remodeling, and membrane/lipid homeostasis with GSH-enhanced detoxification to support biomass recovery and cell integrity under Cr. These 30 candidate hub genes, which are directly or indirectly associated with enhanced photosynthesis, antioxidant capacity, osmotic balance, and growth under Cr stress with mitigator treatment, represent promising targets for sweet potato breeding programs and genetic engineering aimed at improving chromium stress tolerance.

To validate the reliability and consistency of the RNA-seq data, the expression levels of eight candidate genes were analyzed by RT-qPCR using the same CK, Cr, GC, and MC samples (Figure 8). The RT-qPCR expression patterns further corroborated the study’s key functional outcomes. Genes involved in metal transport, such as the metal transporter NRAMP5 (contig2833CG0010), were strongly repressed by Cr (~0.24), and this repression was further enhanced by GSH (~0.13), suggesting that GSH assists the plant in actively limiting Cr uptake. Conversely, genes central to defense and signaling were significantly induced by Cr (> 3.50). The expression of the GST enzyme HSP26-A (scaffold37237CG0020), which is crucial for detoxification, was highly induced by Cr (~6.85) but was notably reduced by GSH treatment (~4.94). This confirms that GSH alleviates the internal oxidative stress load, thereby reducing the necessity for massive GST transcription, whereas MC maintained a high level of induction (~7.64). The transcription factor WRKY70 (scaffold8000.4CG0520) was induced by Cr (~3.54), but its expression was further enhanced by GSH (~4.83), aligning with the finding that GSH actively reinforces stress signaling and defense programs. Furthermore, the MT treatment specifically showed the highest induction for the transcription factor ERF113 (scaffold7262CG0010, approximately 9.05) and the receptor GLR1.3 (~6.73), supporting the role of MT in fine-tuning specific hormone and signaling pathways, as previously indicated by the cluster analysis. These RT-qPCR results collectively validate the RNA-seq dataset and provide a specific functional context for the GSH- and MT-mediated detoxification strategies observed. The correlation between the Log2FC values obtained from the RNA-seq analysis and the RT-qPCR assay was calculated (Figure 8i). The results showed a strong positive Pearson correlation coefficient (r) of 0.915 across the eight selected genes, indicating a high degree of concordance between the two independent quantification methods and robust validation of the transcriptome analysis.