Comparative Transcriptome Analysis of the Molecular Mechanism of the Hairy Roots of Brassica campestris L. in Response to Cadmium Stress

Brassica campestris L., a hyperaccumulator of cadmium (Cd), is considered a candidate plant for efficient phytoremediation. The hairy roots of Brassica campestris L are chosen here as a model plant system to investigate the response mechanism of Brassica campestris L. to Cd stress. High-throughput sequencing technology is used to identify genes related to Cd tolerance. A total of 2394 differentially expressed genes (DEGs) are identified by RNA-Seq analysis, among which 1564 genes are up-regulated, and 830 genes are down-regulated. Data from the gene ontology (GO) analysis indicate that DEGs are mainly involved in metabolic processes. Glutathione metabolism, in which glutathione synthetase and glutathione S-transferase are closely related to Cd stress, is identified in the Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis. A Western blot shows that glutathione synthetase and glutathione S-transferase are involved in Cd tolerance. These results provide a preliminary understanding of the Cd tolerance mechanism of Brassica campestris L. and are, hence, of particular importance to the future development of an efficient phytoremediation process based on hairy root cultures, genetic modification, and the subsequent regeneration of the whole plant.


Introduction
Cadmium (Cd) is one of the most toxic and most common heavy metals, causing serious pollution to farmland and active transfer in soil-plant systems [1][2][3][4][5][6][7][8][9][10][11][12]. Cd intake may lead to tissue inflammation, fibrosis, and even various cancers in humans [13][14][15][16][17], necessitating the development of an efficient strategy to eliminate Cd from soil. Phytoremediation is an effective and low-cost technique for removing Cd from contaminated soil [18][19][20]. Currently, many studies in this area have been carried out using whole plants, which have a limited lifespan and must be replaced and re-established after each experiment. It is also difficult to evaluate the effects of various factors such as light, temperature, soil pH, and rhizosphere microorganisms in different batch experiments [21,22]. Hairy roots, which possess all the morphological and physiological characteristics of normal roots and allow indefinite propagation, have been proven as a convenient experimental tool for investigating the interactions between plant cells and metal ions [22][23][24][25]. The hairy roots of Thlaspi caerulescens were chosen to study oxidative stress and the response of the antioxidant defense system under Cd stress [26]. The hairy roots of Trifolium repens was used as a sensitive tool for monitoring and assessing Cd contamination

An Overview of the mRNA of Brassica campestris L. under Cadmium Stress
To investigate the response of the Cd hyperaccumulator Brassica campestris L. to Cd stress, we analyzed and compared mRNA expression profiles of Cd-treated and untreated hairy roots of Brassica campestris L. Six samples were measured using high-throughput sequencing, and the average data volume of one sample was 6.72 GB (Table 1). Clean reads were obtained from raw reads by removing reads with poor contaminants, low mass, or a high N content of unknown bases. The Q20 value of all clean reads of samples exceeded 98%, and the error rate was less than 1% (Table 1), indicating that the quality of the sequencing met the requirements of subsequent analysis. There were no AT or GC separations in the sequencing report, indicating that sequencing was stable ( Figure S1). Overall, if the proportion of low quality (quality < 20) bases was low, the sequencing quality was good ( Figure S2). After obtaining clean reads, the clean reads were aligned to the reference genome sequence by Hierarchical Indexing for Spliced Alignment of Transcripts (HISAT). The average ratio of each sample reached 64.97% (Table S1), and the uniform ratio between samples indicated that the clean reads data between samples were comparable. To conclude, the transcriptome sequencing results were credible. According to the untreated and treated samples, the similarity heat map shows that the correlations among control replicates C1, C2, and C3 were relatively high; the correlations among replicates within the Cd treatment T1, T2, and T3 were comparatively high. This reveals that Cd treatment caused a large difference in the transcriptome of hairy roots of Brassica campestris L. ( Figure 1A). Through the framed genetic tree ( Figure 1B), it is further illustrated that the closer the expression profile, the higher the correlation, which indicates that Cd stress greatly changed the expression profile of the hairy roots of Brassica campestris L. Above all, the results of transcriptome sequencing were reliable. data volume of one sample was 6.72 GB (Table 1). Clean reads were obtained from raw reads by removing reads with poor contaminants, low mass, or a high N content of unknown bases. The Q20 value of all clean reads of samples exceeded 98%, and the error rate was less than 1% (Table 1), indicating that the quality of the sequencing met the requirements of subsequent analysis. There were no AT or GC separations in the sequencing report, indicating that sequencing was stable ( Figure S1). Overall, if the proportion of low quality (quality < 20) bases was low, the sequencing quality was good ( Figure S2). After obtaining clean reads, the clean reads were aligned to the reference genome sequence by Hierarchical Indexing for Spliced Alignment of Transcripts (HISAT). The average ratio of each sample reached 64.97% (Table S1), and the uniform ratio between samples indicated that the clean reads data between samples were comparable. To conclude, the transcriptome sequencing results were credible. According to the untreated and treated samples, the similarity heat map shows that the correlations among control replicates C1, C2, and C3 were relatively high; the correlations among replicates within the Cd treatment T1, T2, and T3 were comparatively high. This reveals that Cd treatment caused a large difference in the transcriptome of hairy roots of Brassica campestris L. ( Figure  1A). Through the framed genetic tree ( Figure 1B), it is further illustrated that the closer the expression profile, the higher the correlation, which indicates that Cd stress greatly changed the expression profile of the hairy roots of Brassica campestris L. Above all, the results of transcriptome sequencing were reliable.

Analysis of Differentially Expressed Genes (DEGs)
A total of 35,385 genes were detected by sequencing C and T samples. Among these genes, 30,873 genes were found in both C and T samples, while 2490 and 2022 genes were observed in C and T, respectively ( Figure S3). There was a significant difference in unigenes between T and C (p < 0.05), of

Analysis of Differentially Expressed Genes (DEGs)
A total of 35,385 genes were detected by sequencing C and T samples. Among these genes, 30,873 genes were found in both C and T samples, while 2490 and 2022 genes were observed in C and T, respectively ( Figure S3). There was a significant difference in unigenes between T and C (p < 0.05), of which 1564 unigenes were up-regulated and 830 were down-regulated after CdCl 2 treatment (Figure 2.).

Gene Ontology (GO) Functional Analysis of DEGs
Following the identification of DEGs, we analyzed their functions using GO; we classified and enriched their GO function. Detailed proportions of the GO annotation for DEGs are shown in Figure  3, indicating that molecular functions, biological processes, and cellular components were well represented. The GO functional classification results show that the DEGs were mainly enriched in 47 GO terms ( Figure 3). It was found that the biological functions of the DEGs in the hairy roots of Brassica campestris L. were enriched in cell stimulation responding to Cd stress. Among the most abundant genes in hairy roots were those enriched in response to stimulus (43.6%), an organic substance (15.9%), an endogenous stimulus (14.0%), a chemical (21.2%), a carbohydrate (4.6%), stress (21.2%), oxygen-containing (6.4%) and hormone (11.8%) ( Table S2).

Gene Ontology (GO) Functional Analysis of DEGs
Following the identification of DEGs, we analyzed their functions using GO; we classified and enriched their GO function. Detailed proportions of the GO annotation for DEGs are shown in Figure 3, indicating that molecular functions, biological processes, and cellular components were well represented. The GO functional classification results show that the DEGs were mainly enriched in 47 GO terms ( Figure 3). It was found that the biological functions of the DEGs in the hairy roots of Brassica campestris L. were enriched in cell stimulation responding to Cd stress. Among the most abundant genes in hairy roots were those enriched in response to stimulus (43.6%), an organic substance (15.9%), an endogenous stimulus (14.0%), a chemical (21.2%), a carbohydrate (4.6%), stress (21.2%), oxygen-containing (6.4%) and hormone (11.8%) (

Pathway Functional Analysis Results
Different genes in the hairy roots coordinate with each other to perform their biological functions, and pathway-based analysis helps to further understand the biological functions of genes. We mapped the DEGs to the reference canonical pathways in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database to further identify the active metabolic pathways involved in the responses to Cd. The pathway classification results are shown in Figure 4. First, we divided the differentially expressed genes into the following six categories: metabolism, organic systems, environmental information processing, genetic information processing, cellular processes, and human diseases. Then, we further divided the six categories into 20 subcategories. We listed the pathway and the top six differential gene enrichments, as shown in Table S3. It can be seen from the table that Cd stress caused changes in the hairy root metabolic pathway, secondary metabolites, plant-pathogen interactions, plant hormone signal transduction, starch and carbohydrate metabolism, and RNA transport in Brassica campestris L. Altogether, there were 433 significant differentially expressed genes involved in the metabolic signaling pathway, accounting for 24.41% of the total number of differentially expressed genes in the signaling pathway (Table S3).

Real-Time PCR Confirmation of the Gene Expression
To verify the reliability of sequencing, eight differentially expressed genes were randomly selected for RT-PCR validation. Real-time PCR results showed that the transcription levels of these eight genes were consistent with the results of transcriptome sequencing ( Figure 5, Table 2). Using this, we can conclude that the results of transcriptome sequencing can represent differences in transcriptome levels in the hairy roots of Brassica campestris L. caused by Cd stress. We found that the Gene ID 103874543 involved in the synthesis of glutathione was up-regulated.

Pathway Functional Analysis Results
Different genes in the hairy roots coordinate with each other to perform their biological functions, and pathway-based analysis helps to further understand the biological functions of genes. We mapped the DEGs to the reference canonical pathways in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database to further identify the active metabolic pathways involved in the responses to Cd. The pathway classification results are shown in Figure 4. First, we divided the differentially expressed genes into the following six categories: metabolism, organic systems, environmental information processing, genetic information processing, cellular processes, and human diseases. Then, we further divided the six categories into 20 subcategories. We listed the pathway and the top six differential gene enrichments, as shown in Table S3. It can be seen from the table that Cd stress caused changes in the hairy root metabolic pathway, secondary metabolites, plant-pathogen interactions, plant hormone signal transduction, starch and carbohydrate metabolism, and RNA transport in Brassica campestris L. Altogether, there were 433 significant differentially expressed genes involved in the metabolic signaling pathway, accounting for 24.41% of the total number of differentially expressed genes in the signaling pathway (Table S3).

Real-Time PCR Confirmation of the Gene Expression
To verify the reliability of sequencing, eight differentially expressed genes were randomly selected for RT-PCR validation. Real-time PCR results showed that the transcription levels of these eight genes were consistent with the results of transcriptome sequencing ( Figure 5, Table 2). Using this, we can conclude that the results of transcriptome sequencing can represent differences in transcriptome levels in the hairy roots of Brassica campestris L. caused by Cd stress. We found that the Gene ID 103874543 involved in the synthesis of glutathione was up-regulated.

Identification of Key Gene/Protein and Expression Validation
The GO and KEGG analysis of transcriptome sequencing showed that the GSH metabolic pathway plays an important role in regulating the Cd enrichment ability and defense ability of rapeseed hairy roots. Therefore, we continue to study the GSH metabolic pathway. Using the KEGG database, we obtained the following signal path map ( Figure S4). By analyzing key genes in the GSH metabolic signaling pathway, significant up-regulation of GSHB and GST in transcriptome sequencing results caught our attention, so we selected validation of GSHB and GST protein expression levels.
Western Blot indicated that the expression of GSHB at 200 µM Cd stress was 1.85-fold higher than that of the control group, and the expression level of GST was 2.57-fold higher than that of the control group. The amount of actin protein was stable across the treatments ( Figure 6). This result further demonstrated that, not only the mRNA, but also the proteins of GSHB and GST were highly upregulated in the hairy roots of Brassica campestris L.

Identification of Key Gene/Protein and Expression Validation
The GO and KEGG analysis of transcriptome sequencing showed that the GSH metabolic pathway plays an important role in regulating the Cd enrichment ability and defense ability of rapeseed hairy roots. Therefore, we continue to study the GSH metabolic pathway. Using the KEGG database, we obtained the following signal path map ( Figure S4). By analyzing key genes in the GSH metabolic signaling pathway, significant up-regulation of GSHB and GST in transcriptome sequencing results caught our attention, so we selected validation of GSHB and GST protein expression levels.
Western Blot indicated that the expression of GSHB at 200 μM Cd stress was 1.85-fold higher than that of the control group, and the expression level of GST was 2.57-fold higher than that of the control group. The amount of actin protein was stable across the treatments ( Figure 6). This result further demonstrated that, not only the mRNA, but also the proteins of GSHB and GST were highly upregulated in the hairy roots of Brassica campestris L.

Discussion
It is particularly important to investigate the Cd tolerance and the related molecular mechanisms of the hairy roots of Brassica campestris L. [47,50,51]. Our previous studies showed that the high Cd accumulation in the hairy roots of Brassica campestris L. reached more than 1000 mg/kg [47]. The growth of hairy roots of Brassica campestris L. was significantly inhibited when it grow in the medium

Discussion
It is particularly important to investigate the Cd tolerance and the related molecular mechanisms of the hairy roots of Brassica campestris L. [47,50,51]. Our previous studies showed that the high Cd accumulation in the hairy roots of Brassica campestris L. reached more than 1000 mg/kg [47]. The growth of hairy roots of Brassica campestris L. was significantly inhibited when it grow in the medium with CdCl 2 of 200 µM for 24 h [47]. Under this condition, the hairy roots indicated obvious browning and decaying. Additionally, a high Malondialdehyde (MDA) content and microscopic observation of the hairy roots indicated that cell membrane damage and apoptosis appeared, which increased our interest to further investigate the response molecular mechanism using RNA-Seq to perform transcriptome sequencing.
Results indicate that the hairy roots of Brassica campestris L. may improve their Cd tolerance and response to the Cd stress as follows. First, the synthesis of some antioxidants (e.g., GSH and AsA) in the hairy roots may be enhanced by up-regulating expression of the related genes (Figure 7). The KEEG analysis and Western Blot test indicated that the GSHB genes were significantly up-regulated and the GSHB content was improved, hence definitely benefiting the corresponding GSH synthesis ( Figure 6 and Figure S4). When considering GSH as one of the main antioxidants in plant cells [52,53], it, therefore, is expected to enhance the scavenging of reactive oxygen species (ROS) [54].
with CdCl2 of 200 μM for 24 h [47]. Under this condition, the hairy roots indicated obvious browning and decaying. Additionally, a high Malondialdehyde (MDA) content and microscopic observation of the hairy roots indicated that cell membrane damage and apoptosis appeared, which increased our interest to further investigate the response molecular mechanism using RNA-Seq to perform transcriptome sequencing.
Results indicate that the hairy roots of Brassica campestris L. may improve their Cd tolerance and response to the Cd stress as follows. First, the synthesis of some antioxidants (e.g., GSH and AsA) in the hairy roots may be enhanced by up-regulating expression of the related genes (Figure 7). The KEEG analysis and Western Blot test indicated that the GSHB genes were significantly up-regulated and the GSHB content was improved, hence definitely benefiting the corresponding GSH synthesis ( Figure 6 and Figure S4). When considering GSH as one of the main antioxidants in plant cells [52,53], it, therefore, is expected to enhance the scavenging of reactive oxygen species (ROS) [54].
Second, detoxification of the hairy roots from Cd stress could be further enhanced by upregulating the expression of the genes involved in other GSH metabolism aside from synthesis catalyzed by GSHB (Figure 7). The reaction of scavenging ROS can be catalyzed by the enzymes including GPX or GST with GSH as the substrate. Described in the Results section, the GST and GPX (103830386, 103862090) genes in the glutathione metabolic pathway were up-regulated and the enhanced GST synthesis also was verified by a Western Blot test (Figure 6 and Figure S4). Conversely, the up-regulated G6PD (103836988) gene in the Cd-treated hairy roots ( Figure S4) benefitted the synthesis of G6PD, which catalyzed glucose 6-phosphate to produce a co-enzyme (i.e., NADPH) of GSR and help to regenerate GSH from GSSG [55]. Previous studies also observed that the expression of genes related to some key enzymes of GSH metabolism, including ROS scavenging (e.g., GPX and GST) and GSH regeneration from GSSG (e.g., GSR), increased in response to different heavy metal stress (e.g., Ni, Cd etc.) [55][56][57]. Third, the hairy roots of Brassica campestris L. may improve their tolerance to Cd by metalbinding and/or metal-chelating reactions by GSH and/or phytochelatins (PCs). As a major non- Second, detoxification of the hairy roots from Cd stress could be further enhanced by up-regulating the expression of the genes involved in other GSH metabolism aside from synthesis catalyzed by GSHB (Figure 7). The reaction of scavenging ROS can be catalyzed by the enzymes including GPX or GST with GSH as the substrate. Described in the Results section, the GST and GPX (103830386, 103862090) genes in the glutathione metabolic pathway were up-regulated and the enhanced GST synthesis also was verified by a Western Blot test (Figure 6 and Figure S4). Conversely, the up-regulated G6PD (103836988) gene in the Cd-treated hairy roots ( Figure S4) benefitted the synthesis of G6PD, which catalyzed glucose 6-phosphate to produce a co-enzyme (i.e., NADPH) of GSR and help to regenerate GSH from GSSG [55]. Previous studies also observed that the expression of genes related to some key enzymes of GSH metabolism, including ROS scavenging (e.g., GPX and GST) and GSH regeneration from GSSG (e.g., GSR), increased in response to different heavy metal stress (e.g., Ni, Cd etc.) [55][56][57].
Third, the hairy roots of Brassica campestris L. may improve their tolerance to Cd by metal-binding and/or metal-chelating reactions by GSH and/or phytochelatins (PCs). As a major non-enzymatic antioxidant, GSH also may be involved in cell defense against toxicants via binding them or their metabolites with catalysis of GST [53,58]. Previous studies have shown that high levels of expression of most GST genes result in higher GST activity, which increases maize root tolerance to abiotic stress [59]. Data from the GO analysis, regarding the GST enzyme gene in this work, indicate that GST mainly plays a role in the toxin metabolic process and the response to metal ions ( Figures S4 and S5). Alternately, PCs may play an important role in the detoxification of heavy metals in this process [60]. Data analyses indicate that the expression of the glutamate-cysteine ligase catalytic subunit (103861360), which participates in cysteine and methionine metabolism following GSH synthesis and synthesis of subsequent chelating agents (e.g., PCs) with GSH as a precursor, were up-regulated in plants under Cd stress (Figure 7). Previous studies also found that a large number of PCs were observed in rice plants growing under heavy metal stress [53,61]. The functional Vicia sativa Phytochelatin Synthase 1 (PCS1) homolog is overexpressed in Arabidopsis thaliana [62]. The expression levels of GSHB, GSR, PCs, and other chelating peptide synthesis-related enzymes in hyperaccumulators under heavy metal stress were significantly increased [63][64][65][66][67].
To summarize, all these factors could provide an integrated protection for hairy roots from Cd stress. Together with previous studies, results in this work present a better understanding of the regulation network of the hairy roots of Brassica campestris L. under Cd stress, hence providing important information for further genetic modification and subsequent regeneration of the whole plant with improved accumulation abilities.

Plant Material L and Cadmium Stress Treatment
The hairy roots of Brassica campestris L. used in this study were previously induced by infection of the seedlings with Agrobacterium rhizogenes in our laboratory [45]. The hairy roots were maintained in Murashige and Skoog medium (with vitamins and sugar) medium (Coolaber, Beijing, China) based on previous methods [45]. Briefly, 0.5 g of fresh hairy roots of Brassica campestris L. were placed in a 100 mL pyridoxine flask containing 50 mL Murashige and Skoog medium (with vitamins and sugar) liquid medium (NH 4 [47]. Therefore, this treatment condition was chosen for further response mechanism study in this work. The Cd-treated hairy roots of Brassica campestris L. were collected and rinsed three times with deionized water. The samples were then dried and stored at −80 • C. Untreated and Cd-treated samples were labeled C and T, respectively. There were three replicates for each treatment.

Construction and Sequencing of mRNA Libraries
Total RNA was extracted using a TRIzol reagent and was treated with DNase to remove genomic DNA contamination. The eukaryotic mRNA was enriched with Oligo(dT) magnetic beads, the interrupter reagent was used to break the mRNA into short fragments in a Thermomixer (Eppendorf, Hamburg, Germany) to synthesize a strand of cDNA using the interrupted mRNA as a template and, then, the preparation was performed. A two-strand synthesis reaction system synthesized a double-strand cDNA and then used a kit for purification and recovery, cohesive end repair, the addition of base "A" to the 3' end of the cDNA, and the attachment of a linker. Subsequently, the fragment size was chosen, followed by PCR amplification. A suitable library was inspected with an Agilent 2100 Bioanalyzer and an ABI StepOne Plus Real-Time PCR System. Then, we performed paired-end sequencing on an Illumina Hiseq4000 (BGI, Shenzhen, China) following the manufacturer's protocol. The six RNA libraries consisted of three control libraries and three Cd-treated libraries.

Sequence and Primary Analysis
We used the Illumina paired-end RNA-seq approach to sequence the hairy roots of Brassica campestris L. transcriptomes, each producing 6 Gb of multiple data, resulting in a total of 40.32 Gb of sequences. Prior to assembly, low-quality readings were removed, including reads containing sequencing adaptors, sequencing primers, and nucleotides with a quality score below 20. The original sequence data were submitted to the NCBI Database under registration number PRJNA543954.

RNA-Seq Reads Mapping
We compared the readings of the different samples to the Brassica campestris L. (Field mustard) ID: 229 cabbage canola reference genomic sequence using the HISAT2 software. The alignment process can be divided into the following three parts: (1) aligning of reads to a known transcriptome (optional); (2) aligning of the aligned pairs of reads on the reference genome; (3) unpaired read segments are aligned to the reference genome. Reads aligned to the specified reference genome were called mapped reads, and subsequent information analysis was performed based on mapped reads.

Transcript Abundance Estimation and Differential Expression Testing
The results of the HISAT2 [68] alignment to the reference genome were submitted to htseq-count (v0.6.0) [69] for processing, and the read count of each transcript was obtained. Single gene expression levels were calculated using readings per kilobase read (RPKM), which eliminated the effects of gene length and sequencing levels during the calculation of gene expression, making the samples comparable. Then, we used DEGseq for analysis, using the normalization method of quantiles, fold change ≥2, and p-value < 0.05 as the threshold for determining whether the gene was differentially expressed to obtain DEGs.

Gene Annotation, Classification, and Metabolic Pathway Analysis
To study the functional partitioning of DEGs in the hairy roots of Brassica campestris L. under Cd stress, we used GO and KEGG for further annotation, classification, and metabolic pathway analysis [70][71][72]. First, a gene ontology (GO) enrichment analysis of differentially expressed genes was performed by the clusterProfiler R package in which the gene length deviation was corrected. A GO term with a corrected p-value of less than 0.05 was considered to be significantly enriched by differentially expressed genes. The KEGG pathway was retrieved from the KEGG web server (http://www.genome.jp/kegg/). The clusterProfiler R package was used to test the statistical enrichment analysis of differentially expressed genes in the KEGG pathway.

Quantitative Real-Time PCR Validation
Eight DEGs were randomly selected for RT-PCR validation. The primer sequences and reference genes of these genes are listed in Table S4. Total RNA (0.2 µg) from each root sample was reverse transcribed using a PrimScript ® RT Kit (Takara, Beijing, China) and random primers according to the manufacturer's instructions. Quantitative PCR reactions were performed in a 20 µL reaction volume using a Promega GoTaq ® qPCR Master Mix Real Time PCR kit (Takara, Beijing, China) according to the manufacturer's instructions. The reaction was carried out on anSLAN-90P (Hongshi Medical Technology Co., LTD., Shanghai, China). Each biological replica was technically replicated three times. Relative RNA expression of the selected genes, which is the expression of these genes relative to an internal reference gene, was used as an indicator of the genes' expression in different samples [73][74][75]. The relative expression levels of the selected genes were calculated using the 2-∆∆CT method and the probable ubiquitin-conjugating enzyme E2 21 (BnUBC21) gene was used as a reference gene to correct gene expression [76]. Three replicates were performed for each sample.

Western Blot Analysis
The extraction of total hairy root protein was carried out using a Plant Protein Extraction Kit (CoWin Biosciences, Beijing, China). The extracted proteins were assayed by a Pierce™ BCA Protein Assay Kit (Thermo Fisher Scientific, Massachusetts MA, USA). Standard Western blots were performed. Immunoblotting was carried out by using the following primary antibodies: GSH-S (Agrisera AB, Vännäs, Swedish), GST class-phi (Agrisera AB, Vännäs, Swedish), and β-actin (CoWin Biosciences, Beijing, China). After incubation with the secondary antibodies, the signal was developed by chemiluminescence and autoradiography. Densitometric analysis was performed using ImageJ software (National Institutes of Health, Bethesda, NY, USA).

Statistical Analysis
All the experimental data were obtained with three or more repetitions. The data obtained from the experiment were analyzed using IBM's SPSS 20 software by one-way analysis of variance (ANOVA). Statistical analysis was calculated by Duncan's method (p < 0.05).

Conclusions
During this study, 2364 DEGs were discovered under Cd stress in the hairy roots of Brassica campestris L. based on RNA-Seq analysis. These genes were mainly involved in transcription-related processes, defense, stress responses, and transport processes in the response of Brassica campestris L. to Cd stress. Furthermore, data from Western blot tests indicated that the signaling pathway for glutathione synthesis and metabolism played an important role in the response to Cd stress. These results provide valuable information for enhancing our understanding of the heavy metal tolerance of Brassica campestris L. and the corresponding molecular mechanism. It is expected that the Cd-accumulating ability will be further improved by combining transgenic modification of the hairy roots, such as over-expression of the genes involved in Cd hyperaccumulation with subsequent plant regeneration, hence making the whole regenerated plant of Brassica campestris L. more promising for future application in phytoremediation.
Supplementary Materials: Supplementary materials can be found at http://www.mdpi.com/1422-0067/21/1/180/s1. Figure S1: Distribution of base composition on clean reads. X axis represents base position along reads. Y axis represents base content percentage; Figure S2: Distribution of base quality on clean reads. X axis represents base positions along reads. Y axis represents base quality value; Figure S3: Venn diagram analysis; Figure S4: Glutathione metabolism in Brasssica campestris L; Figure S5: Molecular Function with GST unigene; Table S1: Summary of Genome Mapping;

Conflicts of Interest:
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.