Metabolomics analysis reveals different mechanisms of cadmium response and functions of reduced glutathione in cadmium detoxification in the Chinese cabbage

Chinese cabbage (Brassica rapa) is one of the vegetable crops with the largest acreage and yield in East Asia. With the aggravation of heavy metal pollution in soil, the excessive accumulation of heavy metal has become increasingly severe in Chinese cabbage, among which, cadmium is one of the most representative elements. It is critical to understand the regulatory mechanism of Cadmium (Cd) response and tolerance in Chinese cabbage at the metabolism level. Comparative of functional metabolomics and transcript analyses were conducted in Chinese cabbage grown on half-strength MS solid medium with different Cd concentrations (5, 25, 100 µM Cd). The results revealed that the transcript responses were not consistent with the metabolic responses, and Cd treatment at different concentrations could induce more specific differentially expressed metabolites. Moreover, the reduced glutathione (GSH) metabolism pathway was induced under the higher Cd stresses (25, 100 µM Cd). The following result indicated that the GSH played a role in the reduction of Cd accumulation in Chinese cabbage under the higher Cd stresses but led to Cd accumulation under the low Cd stress. This study offers critical information on handling Cd contamination in Chinese cabbages at the different Cd contamination conditions. Different concentration of Cd stresses resulted different responses in Chinese cabbage, and GSH accumulate Cd in the low Cd stress, and exclude Cd in the higher Cd stresses in Chinese cabbage.


Introduction
Due to the improper application of pesticides and chemical fertilizers in recent years, heavy metal pollutants have been increasing in the soil of farm fields and facility agriculture fields, which greatly threats the safety of agriculture production Yang et al. 2019;Zhang et al. 2018;Singh and Prasad 2017). Cadmium (Cd), categorized as one of the class I cancerogens by the World Health Organization (WHO), is a primary element of heavy metal pollutants in soil. Cd mostly exists in the ionic state; most cadmium compounds are easily dissolved in water and absorbed by vegetables and crops, severely increasing the risk of human health through the food chain (Yang et al. 2016).
In plants, some materials and ion transporters that can transport Cd 2+ into the vacuole and out of the cell membrane also function as key Cd detoxification pathways, such as ATP binding cassette (ABC) transporters, heavy metal ATPase (HMA). For example, (heavy metal transporter) HMT1 and HMA3 can isolate Cd into the vacuole in yeast and plants (Huang et al. 2012;Miyadate et al. 2011), and AtPDR8 excludes Cd 2+ outside of the root cell membrane in Arabidopsis (Kim et al. 2007). In Arabidopsis, Nramp1 (natural resistance-associated macrophage protein 1) was identified to mediate the absorbance of Cd and Mn (Cailliatte et al. 2010), and AtHMA2 and AtHMA4 are involved in transporting Cd into xylem (Mills et al. 2005;Wong and Cobbett 2009). In rice, OsNramp1 localizes in plasma membrane Extended author information available on the last page of the article and participates in cellular Cd 2+ transport . OsNramp5 plays an imperative role in Cd and Mn transport in the root of rice, and knockout of OsNramp5 led to the decrease of Cd and Mn accumulation in rice (Sasaki et al. 2012;Yang et al. 2014). OsHMA2 were referred to the loading of Cd into the xylem (Yamaji et al. 2013).
In addition to Cd transporters, plants respond to Cd stress by chelating Cd 2+ with Cd chelating metabolites, thereby reducing the concentration of Cd 2+ and relieving the toxicity of Cd. Reduced glutathione (GSH) is an imperative active material to chelate free radical, which widely exists in plants and animals (Griffith et al. 1978). When plants were under Cd stress, the content of GSH increased rapidly and chelated Cd 2+ , thus relieving the damage of Cd stress to the metabolic activities in the cytosol (Marrs, 1996). Under Cd stress, phytochelatin synthase (PCs) can transfer c-Glu-Cys to GSH, and then generate phytochelatin (PC) (Thangavel et al. 2007;Zhao et al. 2010). PC can not only combine Cd 2+ directly, but also reduce the damage of oxidative stress caused by heavy metals (Schutzendubel et al. 2001;Schutzendubel and Polle 2002). Moreover, GSH and PC facilitate most of ABC transporters in the transport of Cd (Bovet et al. 2005;Rauser 1995;Yamaguchi et al. 2020).
In addition, in plants, metallothionein and organic acid can also bind with Cd 2+ , thereby alleviating the toxicity of Cd to plants (Liu et al. 2020;Choppala et al. 2014;Gallego et al. 2012). In recent years, studies have shown that, under Cd stress, plants generate low molecular weight organic acids, such as oxalic acid and malic acid. These organic acids can be exported from the root, and combine with Cd 2+ , thereby blocking Cd 2+ from the plants, and reducing Cd 2+ absorbance in the root (Montiel-Rozas et al. 2016). Under Cd stress, the root of tomato can secret oxalic acid, thereby keeping Cd 2+ outside, and enhance the Cd tolerance of the tomato (Zhu et al. 2011). Application of citric acid and malic acid to rice roots is able to alleviate the Cd content of the leaves significantly (Sebastian and Prasad 2018). Besides chelating Cd 2+ directly, the metabolites produced by the plants are used to deal with other stresses caused by the heavy metal stresses, such as osmotic and oxidative stress (Arbona et al. 2013;Hernández et al. 2015).
Chinese cabbage (Brassica rapa) is the representative and most popular leafy vegetable in East Asia. Leafy vegetables in the Brassicaceae family have a high ability of Cd accumulation in the leaf parts (Kuboi et al. 1986). With the aggravation of heavy metal pollution in vegetable fields, Cd accumulation in Chinese cabbage is increasingly severe. However, the mechanism associated with Cd accumulation and tolerance is still poorly understood.
Here, we carried out the metabolic and transcript profiling of Chinese cabbage under Cd stresses with different levels of Cd (low (5 μM), medium (25 μM) and high (100 μM)). The analyses of metabolome and transcriptome of Chinese cabbage under different levels of Cd stresses have important theoretical value for us to reveal the mechanism of Cd response and Cd accumulation in Chinese cabbage, which is helpful to cultivate new varieties of low-Cd-accumulation Chinese cabbage.

Plant treatments
Guangdongzao, a cultivar of Chinese cabbage (Brassica rapa), was used in this research, which was widely used in the former traditional Chinese cabbage breeding work, and saved in the Vegetable Research Institute, Shandong Academy of Agricultural Sciences.
The seeds were germinated and grown in half-strength MS solid media with 0, 5, 25 and 100 μM CdCl 2 , respectively, in a controlled environment at 22/20 °C in a 14-h-light/10-h-dark photoperiod for 16 days.
For the measurement of the Cd 2+ concentration in the seedlings, 16-d-old seedlings grown on half-strength MS solid media with different Cd treatments mentioned above were collected and washed as previously described (Gong et al. 2003) with minor modifications. Three biological replicates of each Cd-treated group were collected.
For the measurement of the Cd 2+ concentration in the shoot without or with GSH, the Chinese cabbage seeds were germinated on the soil for 6 days. The seedlings with the same grow condition were transferred to the pots, and then bottom flooded once with 0.4 L per pot (pot size 0.4 L) of 0, 5, 25, and 100 μM CdCl 2 , respectively. After 6 days, 10 milliliter GSH (100 μM) were injected into the soil surrounding the roots, and the seedlings continued growing for 6 days before sampling. Three duplicate samples were taken for each treatment.
The Chinese cabbage cultivar named Beijingxin3 was chose to carry out the measurement of Cd 2+ content in the mature shoots, which is one of the most popular Chinese cabbage cultivar in the market. In this part, two GSH application methods were selected to study the effect of GSH on Cd accumulation in Chinese cabbage with three biological repeats of each experiment. Firstly, GSH was applied by root irrigation. The seeds were sowed at the selected vegetable field, which Cd accumulation was tested. The seedlings were grown in the field for 10 days, and then the soil surround the root was injected with 200 ml water with 0 or 100 μM GSH followed by normal watering until sampling (60 days). Secondly, for foliar GSH application, the seedlings were grown in the field for 30 days, and then 100 ml water with 0 or 100 μM GSH were sprayed in the leaves of the lines, and the Chinese cabbage continued growing for 35 days before three biological samples were collected.
The Cd 2+ contents were measured with an ICP-MC (Thermo iCAP Q, Thermo Fisher Scientific, Waltham, MA, USA). The Cd accumulation of the plant tissues and soil were calculated as the amount of Cd comparing with the fresh weight.

Metabolome with LC-MS/MS and data analysis
Tissues (0.5 g) were grounded with liquid nitrogen and the slurry was resuspended with prechilled 80% methanol and 0.1% formic acid completely. The samples were incubated on ice for 5 min and centrifuged at 15,000 g at 4 °C for 5 min. A part of the supernatant was diluted with LC-MS grade water so that the final concentration of methanol in the solution reached 60%. The samples were subsequently transferred to a fresh tube through a 0.22 μM filter and then centrifuged at 15,000 g at 4 °C for 10 min. Finally, the filtered solution (100 μl) was mixed with pre-cooled methanol (400 μl) and injected into the LC-MS/MS system. LC-MS/MS analyses were performed using a Vanquish UHPLC system coupled with an Orbitrap Q Exactive series mass spectrometer (Thermo Fisher Scientific, Waltham, MA, USA). Samples were injected into a Hyperil Gold column (100 × 2.1 mm, 1.9 μM) using a 16-min linear gradient at a flow rate of 0.2 ml/min. The eluents for the positive polarity mode were solution A (0.1% FA in Water) and solution B (Methanol). The eluents for the negative polarity mode were solution C (5 mM ammonium acetate, pH 9.0) and solution B (Methanol). The solvent gradient was set as follows: 2% B, 1.5 min; 2-100% B, 12.0 min; 100% B, 14.0 min; 100-2% B, 14.1 min; 2% B, 16 min. Q Exactive mass spectrometer was processed in a positive/negative polarity mode with a spray voltage of 3.2 kV, a capillary temperature of 320 °C, a sheath gas flow rate of 35 arb, and an aux gas flow rate of 10 arb.
The raw data generated by UHPLC-MS/MS were processed using the compound discoverer 3.0 software (v. 3.0, Thermo Fisher Scientific, Bremen, Germany) to perform peak alignment, peak picking, and quantitation for each metabolite. The parameters were set as follows: retention time tolerance, 0.2 min; actual mass tolerance, 5 ppm; signal intensity tolerance, 30%; signal/noise ratio, 3; and minimum intensity, 100,000. After that, the peak intensities were normalized to the total spectral intensity. The normalized data predict the molecular formula based on additive ions, molecular ion peaks and fragment ions. The peaks were matched with the mzCloud (https:// www. mzclo ud. org/) and ChemSpider (http:// www. chems pider. com/) database to obtain the accurate qualitative and relative quantitative results. Statistical analysis was performed using the statistical software R (R version R-3.4.3), Python (Python 2.7.6 version) and CentOS (CentOS release 6.6). When data were distributed and disordered, normal transformations were attempted using the area normalization method.
These metabolites were annotated using the KEGG database (http:// www. genome. jp/ kegg/), HMDB database (http:// www. hmdb. ca/) and Lipid map database (http:// www. lipid maps. org/). Principal components analysis (PCA) and Partial least squares discriminant analysis (PLS-DA) were conducted using the metaX software (Wen et al. 2017). We applied univariate analysis (t-test) to calculate the statistical significance (P value). The metabolites with VIP > 1, P value < 0.05 and fold change (FC) ≥ 1.2 or FC ≤ 0.833 were considered to be differentially expressed metabolites. Volcano plots were used to filter metabolites of interest based on log 2 (FC) and − log 10 (P value) of the metabolites. For clustering analysis, the data were normalized using z-scores of the intensity areas of the differentially expressed metabolites, and plotted using the P heat map package. The correlation between the differentially expressed metabolites was analyzed by cor () in the R language (method = pearson). A statistically significant correlation between the differentially expressed metabolites was calculated by cor. m test () in the R language. P value < 0.05 was considered as statistically significant, and correlation plots were plotted by the corrplot package in the R language. The functions of these metabolites and metabolic pathways were analyzed using the KEGG database. The significant metabolic enrichment pathways of differentially expressed metabolites were investigated, and the P value of the metabolic pathways < 0.05 was selected as the differentially regulated pathways.

Transcriptome sequencing and data analysis
Three biological replicates selected in half-strength MS solid media without or with Cd treatments were used to perform RNA-seq analysis. The total RNA was isolated by Trizol, and mRNA was purified from total RNA using poly-T oligoattached magnetic beads. First-strand cDNA was synthesized using a random hexamer primer and M-MuLV Reverse Transcriptase (RNase H-). Second strand cDNA synthesis was subsequently performed using DNA polymerase I and RNase H. The library fragments were purified with the AMPure XP system (Beckman Coulter, Beverly, USA). Then PCR was performed with Phusion High-Fidelity DNA polymerase, universal PCR primers and index (X) Primer. Finally, PCR products were purified (AMPure XP system), and the library quality was evaluated on the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA).
The clustering of the samples was conducted on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumia, San Diego, USA) according to the manufacturer's instructions. After cluster generation, the cDNA libraries were sequenced on an Illumina Novaseq platform to generate 150 bp paired-end reads. For novel transcript prediction, each sample's mapped reads were assembled by StringTie (v1.3.3b) (Pertea et al. 2015). StringTie used a novel network flow algorithm and an optional de novo assembly step to assemble and quantitate full-length transcripts, representing multiple splice variants for each gene locus. Novogene carried out total RNA isolation, library construction, sequencing, and raw data analysis.

Correlation analysis between DEMs and DEGs
Three independent biological replicates from the six were randomly selected in the DEMs and co-analyzed with DEGs from the independent biological replicates. Fifty DEMs and one hundred DEGs whose P value ranged from the lowest value among all the DEMs and DEGs were analyzed by the Pearson correlation coefficient (Adler and Parmryd 2010).

Statistical analysis
All experiments were conducted in triplicates. Data were analyzed using variance (ANOVA) analysis with Duncan's multiple range test (DMRT). Standard errors were calculated for all mean values, and differences were considered significant at the P ≤ 0.05 level.

Phenotypic variation of Chinese cabbage under different levels of Cd stresses
As shown in Fig. 1A and B, the seedlings treated with the medium-Cd stress showed a significantly worse growth condition and a shorter root length compared with those of the seedlings treated with the low-Cd-stress, and the seedlings grown under high-Cd-stress were also significantly worse than the seedlings under the medium-Cd-stress. Under the high Cd concentration, the Cd accumulation in the Chinese cabbage was greatly elevated (Fig. 1C). These results suggested that stresses with different Cd concentrations could cause different damages for Chinese cabbage and induce different reactions in response to the stresses.

Differentially expressed metabolites (DEMs) identified from metabolomics analysis
To detect the difference in metabolites and pathways in response to different Cd stresses, the metabolic profiles in the 16-d-old seedlings without or with the Cd treatments were analyzed by LC-MS in the positive ion detection model (ESI+) and the negative ion detection model (ESI −). A total of 2275 and 903 metabolites were detected in all the samples in ESI + and ESI −, respectively (Supplemental Table S1). Principal component analysis (PCA) was performed to confirm the quality of the data from the different samples. As shown in Fig. 2A and B, the control sample clustered together and all six independent biological repeats from each group surround the quality control samples. The partial least-squares discriminant analysis (OPLS-DA) was also used to detect the differences between the Cd treatment groups and the control group ( Fig. 2C-I). The DEMs comparing Cd-treated plants against the control were selected according to variable importance in the projection (VIP) based on the selection criteria of OPLS-DA > 1, fold change (FC) > 1.2 (or FC < 0.833), and P value < 0.05 (t-test) (Supplemental Table S2-8).

Profiling of DEMs expressed under all three Cd stresses in Chinese cabbage
To better and concisely understand the difference, the metabolic profile of 5 μM Cd-treated plants compared to the control is hereafter referred to as comparison 1 or C1; likewise, 25 μM Cd-treated plants compared to the control is C2; and 100 μM Cd-treated plants compared to the control is C3. As shown in Fig. 3A, there were 157 DEMs (101 were up-regulated and 56 were down-regulated) in C1, 177 DEMs (108 were up-regulated and 69 were down-regulated) in C2, and 297 DEMs (141 were up-regulated and 156 were downregulated) in C3. Seven were common to both C1 and C2, 29 were common to both C1 and C3, 19 were common to both C2 and C3, and 18 were common to all the comparison groups (Fig. 3C). The total number of the DEMs is 509, and the majority of these DEMs had different fold changes in expression (Fig. 3B). Detailed information regarding these metabolites was summarized in Supplemental Table S9. Among the 509 DEMs, only 3.5% (18) were found to be common in all three comparisons including methanediylidenebis, FzM1, virstatin, convicine, doramapimod, dictyoquinazol A, regorafenib and etc. (Supplemental Table S8).

DEMs specifically regulated in three different Cd stresses in Chinese cabbage
In addition to the 18 common DEMs (3.5%) in all three comparisons, 7 DEMs (1.4%) were involved in both C1 and C2; 29 DEMs (5.7%) were involved in both C1 and C3; and 49 DEMs (9.6%) were involved in both C2 and C3. By the intersection analysis (Fig. 3C), a significant number of DEMs were found to be specifically regulated in C1, C2 and C3, respectively. In C1, 20.2% of the selected DEMs (103) were found to be specifically regulated. In C2, the same percentage of the DEMs (103) was found to be specifically regulated. In C3, 39.3% of the DEMs (200) were specifically involved. Compared to the DEMs in the cross sections, there were many more DEMs specifically induced under the low-Cd-stress, the medium-Cd-stress, and the high-Cd-stress treatments.

Metabolic pathways corresponding to different Cd stresses
To better understand the responses of Chinese cabbage to low-Cd-stress, medium-Cd-stress, and high-Cd-stress, the DEMs were mapped to the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. As shown in Table 1, with low-Cd-stress, phosphonate and phosphinate metabolism and indole alkaloid biosynthesis pathways were differentially induced. When treated with medium-Cd-stress and the high-Cd-stress, only glutathione metabolism pathway was differentially regulated. However, more DEMs were involved in glutathione metabolism pathway in high-Cd-stress (homotrypanothione, glutathionylspermine, glutathione disulfide, glutathionylaminopropylcad averine) than in medium-Cdstress (glutathionylspermine, glutathione disulfide), and the P value of the pathway was much smaller under the  high-Cd-stress than the medium-Cd-stress (Table 1), indicating that glutathione metabolism pathway was more significantly regulated in response to the high-Cd-stress compared with the medium-Cd-stress.

Transcriptome analysis in the Chinese cabbage under the three different Cd stresses
The different growth situation and Cd accumulation of the Chinese cabbage under the treatment of the three Cd concentrition could also be associated with gene expression. RNA-seq analysis was conducted on the three comparisons (C1, C2 and C3) to verify the DEGs' expression pattern. Differentially expressed genes (DEGs) with statistically significant change (a fold change > or < 1, DESeq adjusted P value < 0.05) (Love et al. 2014) in C1, C2, and C3 were selected. In the three different treatments, the largest number of DEGs (2288 up-regulated and 1575 down-regulated) was identified in C3 (100 μM Cd vs control), and the least number of DEGs (350 up-regulated and 539 down-regulated) was identified in C1 (5 μM Cd vs control) (Fig. 4A). As shown in Fig. 4B, 150 DEGs are common to the three comparisons. However, more DEGs specifically exist in different comparisons, and the number of DEGs enhance as the increase of the Cd concentration treatment.
To validate the RNA-seq data, 7 DEGs in the intersection part of the three comparison groups, which shows a log 2 (Fold Change) ≥ 2, were selected for real-time quantitative PCR (qPCR) identification. As shown in Fig. 4C, most of these DEGs share similar expression patterns by the two analyses, indicating that RNA-seq could be used for gene expression profiling of Chinese cabbage in response to the different levels of Cd stress. Among these DEGs, BraA05g004360.3C and BraA04g002030.3C are members of glycoside hydrolase family encoding proteins similar to betaglucosidase, which had a peak value at a certain concentration of Cd stress. Another two DEGs, BraA08g023600.3C and BraA03g002980.3C encoding a gibberellin 2-oxidase and a methionine sulfoxide reductase 3, respectively, were enriched to the same GO term (GO: 0,055,114, oxidation-reduction process). BraA08g015710.3C is a member of SAUR family protein, relating to plant hormone signal transduction, while BraA04g026280.3C and BraA01g019180.3C are unknown genes.

Co-analysis of metabolomics and transcriptome
Correlation analysis was performed between transcriptomics and metabolomics to reveal the consistency between DEGs and DEMs. One hundred DEGs and fifty DEMs were randomly selected according to the P value (from lowest to highest scores) and compared. As shown in Supplementary  Fig.S1, under lower Cd stresses (C1 and C2), the correlations between DEGs and DEMs are over eighty percentage; under higher Cd stress (C3), the correlations is just about fifty percentage. It is suggested that the response of Chinese cabbage to Cd stress was deeply interrupted and more complicated at the transcriptional and metabolic level, and the DEGs did not directly associate with the DEMs under severe Cd stress. According to the KEGG database, the DEMs and DEGs were compared and assigned to explore the fundamental functional cadmium-response pathways further. As shown in Table 1 and 2, three KEGG pathways were identified in the four comparison groups at the metabolic level, and sixteen KEGG pathways were identified at the transcriptional level. We also performed a correlation analysis of KEGG pathways at both metabolic and transcriptional levels ( Supplementary Fig.S2). The results showed that only two pathways were same enriched in the metabolic and transcriptional levels.

GSH accumulated Cd in the Chinese cabbage with the low-Cd-stress and reduced Cd content with the medium-and high-Cd-stress
As the metabolites are the ultimate products of the biological process, and the direct substrates of stress responses, and some of the metabolites could be used for improving the growth state on the adversity stress, we focused on the result of the metabolic analysis. As shown in Fig. 5A, under the medium-Cd stress, two metabolites, glutathione disulfide (GSSG) and glutathionylspermine, were significantly reduced compared to the control. Under the high-Cd-stress, more two metabolites: glutathionylaminopropylcadaverine and homotrypanothione were also significantly declined. Considering all of these four metabolites were mainly generated by GSH in the GSH metabolism pathway, GSH might play a critical role in Cd responses in Chinese cabbage under Cd stresses at higher concentration.
To further identify the role of GSH in Cd responses in Chinese cabbage, we germinate the "Guangdongzao" seeds for 6 days, and then transfered the seedlings into the pot with soils from the vegetable farmland (36º71′43.83″ N, 117º09′08.62″ E) (Table 3), which bottom flooded with 5, 25 and 100 μM Cd solution. The exogenous GSH (100 μM, 10 ml) were added in the soil surrounding the roots of Chinese cabbage seedlings after 6 days. The shoot was collected for the Cd content measurement after another 6 days. As shown in Fig. 5B, under the low-Cd-stress, GSH facilitated the Cd accumulation in the shoot of Chinese cabbage. However, under the medium-Cd and high-Cd-stress, GSH reduced the Cd accumulation.

GSH elevated the Cd accumulation in the shoot of the Chinese cabbage in the vegetable lands
To test whether GSH can be utilized for alleviating the Cd contamination in the production of Chinese cabbages, we sowed the seeds of the representative Chinese cabbage cultivar (Beijingxin3) in China in the vegetable fields, which the Cd content of the soil was below the Chinese national standard (36º71′43.83″ N, 117º09′08.62″ E) (Table 3) or over the standard (36º78′51.33″ N, 118º84′24.81″ E) ( Table 4). As shown in Fig. 5C and D, when treated with a certain concentration (100 μM) of GSH, the Chinese cabbages grown on the soils with different Cd contents accumulated more Cd than the control. For foliar GSH application, the GSH also can accumulate Cd in the leaves of Chinese cabbage, and reduce the Cd contamination in the soil (Fig. 5E). This indicated that, in the tested vegetable field, GSH enhanced the Cd accumulation in the Chinese cabbages.

Discussion
Stresses with different Cd concentrations could cause different damages for Chinese cabbage and induce different reactions in response to the stresses. To clarify the differences in metabolite and transcript pathways in response to different Cd stresses, comparative of functional metabolite and transcript analyses were carried out in Chinese cabbage treated with different Cd concentrations.
Among the 509 DEMs detected in the metabolite analysis, only 3.5% were found to be common in all three comparisons, which indicated that these common DEMs were likely induced directly by the Cd toxicity, and might have little relation with the concentration of Cd. And the results also suggested that notable metabolic differences exist in Chinese cabbage under different Cd stress levels.
Similarly, as the concentration of the cadmium stress increased, the expression of the genes in Chinese cabbage was affected severely, and the different Cd stresses of Chinese cabbage could share specific cadmium response mechanisms at the transcriptional level. Correlation analysis of KEGG pathways was carried out at both metabolic and transcriptional levels, we found that only two pathways were same enriched, suggesting DEMs and DEGs' enrichments from the KEGG database were not consistent in the Chinese cabbage with cadmium stress. Considering metabolites, as the ultimate products, play important roles in plant growth and development, we focused on the metabolite analysis.
KEGG analysis on metabolomics showed that only glutathione metabolism pathway was differentially regulated when treated with medium-and high-Cd-stress, which suggested that the Cd stresses at different concentrations caused different responses in Chinese cabbage, and glutathione metabolism pathway was critical in response to the higher Cd stresses.
In the following test, we found that GSH reduced Cd accumulation in the Chinese cabbage under the medium and high Cd stresses, however, GSH elevated the Cd accumulation under the low Cd stress. It suggested that GSH could play different roles in Cd accumulation of Chinese cabbage under the different Cd concentration. GSH is a wellknown Cd chelation compound, which can bind with Cd 2+ and other heavy metal ions (Jozefczak et al. 2014;Rizwan et al. 2018). In plants, GSH facilitates Cd transporters in transporting Cd 2+ , which can be used to isolate or remove Cd 2+ to prevent the metabolism in the cell (Brunetti et al., 2015;Rauser, 1995).
However, in our field test, the Cd accumulation in the Chinese cabbages was enhanced after GSH treatment. This is different from reports about the function of GSH in the 1 3

BraA01g019180.3C
∇ ∇ prevention of Cd contamination in rice and oilseed rape (Cai et al. 2010;Cao et al. 2015;Huang et al. 2019). Considering the mild Cd contamination could be the popular Cd contamination occurred in the vegetable field (Kumar et al. 2007;Ma et al. 2018), our results indicated that GSH could not be used as Cd scavenger in the most vegetable field for the production of Chinese cabbages.

The proposed model
Considering GSH facilitates Cd transporters in Cd transporting, the difference of Cd accumulation could be due to the different kinds of Cd transporters. As shown in Fig. 5, we speculate that, under the low-Cd-stress, Cd 2+ was isolated by the Cd transporters in the vacuole, which might be the main Cd tolerance pathway. With the GSH, more Cd 2+ were transported to the vacuole from cytosol, and more Cd 2+ will enter into the cytosol from the outside of the root by the ion gradient, which is consistent with the observation in rice (Miyadate et al. 2011). This tolerance mechanism will relieve the damage of Cd in the cytosol. However, more Cd will be accumulated in the whole plant. When Chinese cabbage was treated under medium-and high-Cd-stress, Cd exclusion played the main role in the resistance of Cd toxicity; with GSH, more Cd 2+ were transported outside of the cell by the Cd transporters in the root cell membrane, which has been observed in Arabidopsis (Kim et al. 2007). Moreover, GSH significantly affect iron nutrition, and the real function of most "Cd transporters" is the transporters of some iron ions. Thus, in fact, Cd 2+ is the competitor of the iron ions. According to our result, GSH facilitate  Green boxes indicate metabolites with decreased abundance, while blank boxes imply no significant changes, and the number of the inside boxes represent the fold change of the metabolite in C2 or C3. In each box, the left compartment represents metabolite alteration in C2 and the right compartment indicates the corresponding alteration in C3. B Cd content in shoot in 18-day-old Chinese cabbage (Guangdongzao) seedlings. GSH (100 µM, 10 ml) was applied until they were 12 days old. C, D Cd content of the 65-day-old mature tissues of the Chinese cabbage (Bei-jingxin3) in the soil which contain 0.385 ± 0.020 mg/kg Cd 2+ (C) and 0.089 ± 0.010 mg/kg (D). The seedlings were grown in the field for 10 days, and then the soil surround the root was injected with 200 ml water without or with 100 μM GSH followed by normal watering until sampling. Dotted and segment lines indicate Chinese national standard for cadmium content in agricultural land and vegetables respectively. P values from Student's t test were determined for plants treated with Cd and GSH compared with plants treated with GSH: *P < 0.05; **P 0.01. E Cd content of the Chinese cabbage planted in the Cd contaminated soil (0.385 ± 0.020 mg/kg Cd 2+ ) and the soil surround the Chinese cabbage after the treatment of water and GSH. The seedlings were grown in the field for 30 days, and then 100 ml water without or with 100 μM GSH were sprayed in the leaves of the lines, and the Chinese cabbage continued growing for 35 days before sampling. P values from Student's t test were determined for plants treated with GSH compared with plants treated with water: *P < 0.05; **P 0.01. FW means Fresh weight. Error bars indicates ± SD from three or five (Cd content of soil) independent experiments ◂ Data availability The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request. RNA-Seq sequence data were deposited in the NCBI Sequence Read Archive (SRA, http:// www. ncbi. nlm. nih. gov/ Traces/ sra) under accession number SAMN22567713, SAMN22567714, SA M N 2 2 5 6 7 7 1 5 , SA M N 2 2 5 6 7 7 1 6 , SA M N 2 2 5 6 7 7 1 7 , SA M N 2 2 5 6 7 7 1 8 , SA M N 2 2 5 6 7 7 1 9 , SA M N 2 2 5 6 7 7 2 0 , SA M N 2 2 5 6 7 7 2 1 , SA M N 2 2 5 6 7 7 2 2 , SA M N 2 2 5 6 7 7 2 3 , SAMN22567724 for C_1, C_2, C_3, VP25_1, VP25_2, VP25_3, VP50_1, VP50_2, VP50_3, VP100_1, VP100_2, and VP100_3, respectively.
Code availability Not applicable.