- Research article
- Open Access
Comparative transcriptome analysis reveals key cadmium transport-related genes in roots of two pak choi (Brassica rapa L. ssp. chinensis) cultivars
BMC Genomics volume 18, Article number: 587 (2017)
Cadmium translocation from roots to shoots is a complex biological process that is controlled by gene regulatory networks. Pak choi exhibits wide cultivar variations in Cd accumulation. However, the molecular mechanism involved in cadmium translocation and accumulation is still unclear. To isolate differentially expressed genes (DEGs) involved in transporter-mediated regulatory mechanisms of Cd translocation in two contrasting pak choi cultivars, Baiyewuyueman (B, high Cd accumulator) and Kuishan’aijiaoheiye (K, low Cd accumulator), eight cDNA libraries from the roots of two cultivars were constructed and sequenced by RNA-sequencing.
A total of 244,190 unigenes were obtained. Of them, 6827 DEGs, including BCd10 vs. BCd0 (690), KCd10 vs. KCd0 (2733), KCd0 vs. BCd0 (2919), and KCd10 vs. BCd10 (3455), were identified. Regulatory roles of these DEGs were annotated and clarified through GO and KEEG enrichment analysis. Interestingly, 135 DEGs encoding ion transport (i.e. ZIPs, P1B-type ATPase and MTPs) related proteins were identified. The expression patterns of ten critical genes were validated using RT-qPCR analysis. Furthermore, a putative model of cadmium translocation regulatory network in pak choi was proposed.
High Cd cultivar (Baiyewuyueman) showed higher expression levels in plasma membrane-localized transport genes (i.e., ZIP2, ZIP3, IRT1, HMA2 and HMA4) and tonoplast-localized transport genes (i.e., CAX4, HMA3, MRP7, MTP3 and COPT5) than low Cd cultivar (Kuishan’aijiaoheiye). These genes, therefore, might be involved in root-to-shoot Cd translocation in pak choi.
Soil contamination with cadmium (Cd) has long been a major ecological concern worldwide because of the potential threat to all organisms. Cd is readily transferred from contaminated soil to plants and accumulates in different organs [1,2,3]. The accumulation of Cd in plants may lead to severe toxicity in both plants and animals/humans . To prevent this negative effect, two strategies have been proposed: (i) reducing Cd concentration in contaminated soils using phytoremediation and chemical remediation [4, 5]; and (ii) reducing Cd concentration in edible parts of plants by breeding low-Cd accumulation cultivars. For these purposes, understanding the molecular mechanisms underlying Cd uptake, translocation and accumulation in plants is critically important.
Several processes are involved in the distribution and translocation of Cd in plant roots [6,7,8], including (i) transporter-mediated Cd translocation via symplasmic or apoplastic pathway; (ii) deposition of Cd in the cell wall; and (iii) sequestration of Cd-chelates in the vacuole. Most of these processes are mediated by transport-related genes. In rice, Satoh-Nagasawa et al.  report that loss of OsHMA2 function in insertion mutant results in decreased leaf and grain Cd concentrations. OsHMA3 plays a role in the sequestration of Cd into vacuoles in root cells [10, 11], and OsHMA3 expression in the Cd-sensitive Δycf1 mutant can cause a wild-type tolerance . Furthermore, the member of Nramp family in rice, OsNramp1  and OsNramp5 [13,14,15] were reported to be Cd transporter involved in the root uptake of Cd from the medium. These findings suggested that a number of metal-regulated transporter proteins could participate in cellular Cd uptake and translocation within plants.
High-throughput cDNA sequencing (RNA-seq) is a recently developed approach to transcriptome profiling that has been used to discover and characterize genes, analyze functional gene variation, and identify and quantify rare transcripts [16, 17]. Using this approach, many Cd stress-related genes have been identified in Arabidopsis , ramie , radish  and Viola yedoensis . In addition, comparative transcriptomic analyses have been used for revealing the mechanisms of Cd accumulation in plants [21,22,23]. These studies have greatly improved our understanding of the molecular regulation mechanisms underlying Cd tolerance and accumulation.
Pak choi (Brassica rapa L. ssp. chinensis) is an important leafy vegetable crop. Variability among pak choi cultivars in Cd concentration has been reported [24, 25]. However, the molecular mechanism of Cd accumulation in pak choi is not well understood . Although a fraction of DEGs related to Cd stress has been identified in pak choi , it might not broadly and deeply reveal the mechanism of Cd uptake and translocation.
Here, a comparative transcriptome analysis was performed using RNA-seq technology. The aims were: (i) to obtain assembled unigenes from pak choi root transcriptome; (ii) to investigate the gene expression patterns in the roots of two distinct pak choi genotypes under long-term Cd exposure; and (iii) to identify DEGs involved in Cd-transport regulatory networks.
Illumina sequencing and transcriptome assembly
A total of 119.6 × 106, 114.5 × 106, 105.8 × 106 and 118.3 × 106 raw reads were generated from two replicate libraries for BCd0, BCd10, KCd0 and KCd10, respectively (Table 1). After filtering the reads with low quality and adapters, the percentages of clean reads in all eight transcriptomes were all above 95.99%, and the proportion of Q20 bases for the clean reads was above 95.93% for each library (Table 1). Using the Trinity program, a total of 334,411 putative transcripts were obtained, with an average length of 718 bp and an N50 of 1321 bp, and transcripts with lengths of more than 500 bp accounted for about 37.78% of all transcripts (Fig. 1a and c). After comparing the different transcripts representing one unigene, the longest transcript for each locus was selected as the unigene. In all, 244,190 unigenes were obtained as reference transcripts of pak choi. The mean length was 512 bp, and unigenes with lengths of more than 500 bp accounted for about 23.58% of all unigenes (Fig. 1b and c).
Functional annotation and classification of non-redundant unigenes
To investigate potential functions of the assembled unigenes, sequence similarity searches of 244,190 unigenes were performed in the public databases. In total, 142,631 unigenes were annotated representing 41.59% of the assembled unigenes. In the Nr, Pfam, GO, KOG and KEGG databases, 108,001, 99,209, 100,852, 61,894 and 53,311 unigenes were aligned, respectively (Table 2).
GO assignments were used to classify the functions of all predicted unigenes based on the annotations from the Nr and Pfam databases. In total, 100,852 unigene sequences (41.3%) were categorized into 56 functional groups consisting of 25 biological process, 21 cellular components and 10 molecular function subcategories (Fig. 2a).
The sequence similarity search was performed against the KOG databases to obtain the functional annotations of assembled unigenes. 61,894 unigene sequences (25.34%) with significant homology were assigned to 26 KOG categories (Fig. 2b). Among them, the five largest groups included ‘translation, ribosomal structure and biogenesis’ (9479, 15.31%), ‘posttranslational modification, protein turnover, chaperones’ (9004, 14.55%), ‘general function prediction’ (7987, 12.90%), ‘signal transduction mechanisms’ (6279, 10.14%) and ‘energy production and conversion’ (3890, 6.28%).
To further characterize molecular-level gene functions in terms of biological system networks, the assembled unigenes were mapped against the KEGG protein database. 53,311 unigenes with significant matches were assigned to 132 KEGG pathways, of which the top five were ‘ribosome’, ‘carbon metabolism’, ‘biosynthesis of amino acids’, ‘protein processing in endoplasmic reticulum’ and ‘RNA transport’ (Additional file 1: Table S1). These pathways belonged to five main categories that could be further divided into 19 subcategories (Fig. 3). Among them, translation was the dominant pathway (9048 unigenes, 16.97%), followed by carbohydrate metabolism (5031, 9.44%), overview (4535, 8.51%), folding, sorting and degradation (4372, 8.20%) and transport and catabolism (3778, 7.09%).
Identification and functional annotation of DEGs
A total of 6827 unigenes were found to have significant differential expression in the four comparisons, including BCd10 vs. BCd0 (690 unigenes), KCd10 vs. KCd0 (2733 unigenes), KCd0 vs. BCd0 (2919 unigenes) and KCd10 vs. BCd10 (3455 unigenes) (Additional file 2: Table S2; Fig. 4a). The Venn diagram indicates the distribution of DEGs among the four comparisons (Fig. 4b). Among of these DEGs, up-regulated DEGs accounted for 53.33% of BCd10 vs. BCd0, 93.49% of KCd0 vs. BCd0, and 54.44% of KCd10 vs. BCd10 DEGs, whereas only 3.62% of DEGs were up-regulated in KCd10 vs. KCd0 (Fig. 4a ). There were 2239 DEGs that showed similar expression patterns in the two cultivars, while 2076 DEGs showed opposite expression patterns (Additional file 2: Table S2). Furthermore, 2265 DEGs were specifically regulated by Cd in Kuishan’aijiaoheiye, while only 137 specific regulated DEGs were found in Baiyewuyueman (Additional file 2: Table S2).
GO and pathway enrichment analysis were conducted to gain functional annotations of the DEGs. A total of 70, 67, 93 and 46 terms were significantly enriched for DEGs in BCd10 vs. BCd0, KCd10 vs. KCd0, KCd0 vs. BCd0, and KCd10 vs. BCd10, respectively, and of these, 62.86%, 2.99%, 83.87% and 76.09%, respectively, were up-regulated (Additional file 3: Table S3). In Baiyewuyueman, the predominant enriched GO terms of the Cd-induced up-regulated DEGs were related to cellular process (GO: 0009987), metabolic process (GO: 0008152) and organic substance metabolic process (GO: 0071704) in biological processes. Among the down-regulated DEGs, the mainly enriched GO terms were related to transcription, DNA-templated (GO: 0006351), nucleic acid-templated transcription (GO: 0097659) and RNA biosynthetic process (GO: 0032774) in biological processes. In Kuishan’aijiaoheiye, the majority of annotated DEGs (97.01%) were down-regulated by Cd. The predominant enriched GO terms of these DEGs were involved in biological processes including cellular process (GO: 0009987), single-organism cellular process (GO: 0044763) and cell (GO: 0005623). Compared with Baiyewuyueman, some predominantly enriched GO terms of down-regulated DEGs related to biological processes such as response to oxidative stress (GO: 0006979), cell wall organization (GO: 0071555), glucan metabolic process (GO: 0044042) and obsolete peroxidase reaction (GO: 0006804), were observed in Kuishan’aijiaoheiye under Cd treatment (Additional file 3: Table S3).
Additionally, a total of 44 and 40 significantly enriched pathways were obtained for up and down-regulated DEGs, respectively (Additional file 4: Table S4). As expected, the DEGs annotated in the pathway showed different effects in Baiyewuyueman and Kuishan’aijiaoheiye exposed to Cd. The most significantly enriched pathway annotations for Baiyewuyueman were up-regulated DEGs under Cd stress, such as ‘ribosome’ (74 unigenes), ‘legionellosis’ (10 unigenes), ‘antigen processing and presentation’ (7 unigenes) and ‘carbon fixation in photosynthetic organisms’ (6 unigenes) etc. (Additional file 4: Table S4). However, the most significantly enriched annotations for Kuishan’aijiaoheiye were down-regulated DEGs under Cd stress. The predominantly enriched pathways for the down-regulated DEGs were involved in ‘purine metabolism’ (41 unigenes), ‘spliceosome’ (29 unigenes), ‘rap1 signaling pathway’ (29 unigenes), ‘oxytocin signaling pathway’ (29 unigenes) and ‘leukocyte transendothelial migration’ (29 unigenes) (Additional file 4: Table S4). Furthermore, some Cd stress-related pathways, such as glutathione metabolism (ko00480), ascorbate and aldarate metabolism (ko00053), glycine, serine and threonine metabolism (ko00260) and plant hormone signal transduction (ko04075) are significantly enriched for the up-regulated DEGs in KCd10 vs. BCd10.
DEGs involved in Cd transport
In the four comparisons, 135 DEGs were identified as having high similarity with diverse transporters such as zinc transporters (e.g. ZIP2, ZIP3 and ZIP10), Fe2+ transport protein 1 (IRT1), P1B-ATPase (e.g. HMA2, HMA3, HMA4 and HMA5), ABC transporters, vacuolar cation/proton exchanger 4 (CAX4), oligopeptide transporter (OPT) families and metal tolerance proteins (MTPs) (Fig. 5; Additional file 5: Table S5). Of 135 DEGs, 11 and 4 unigenes were up-regulated by Cd, while 5 and 44 unigenes were down-regulated in Baiyewuyueman and Kuishan’aijiaoheiye respectively. All unigenes belonging to ZIPs, high affinity nitrate transporter 2.6 (NRT2.6) and HMA3 were tended to be up-regulated in both cultivars compared to their respective controls. All unigenes encoding CAX4 and cyclic nucleotide gated channels (CNGCs) were up-regulated in Baiyewuyueman and down-regulated in Kuishan’aijiaoheiye compared to their respective controls. Based on pairwise comparisons, members within the some families such as OPTs, ATP synthases, V-type proton ATPase (V-ATPase), P-type ATPase superfamily (P-ATPase), inorganic pyrophosphatase (PPase), cadmium/zinc/copper-transporting ATPase (HMAs) and metal tolerance proteins (MTPs), had considerably variable in expression pattern in the four groups. Moreover, the same member within the same pairwise comparison also had a different expression pattern; for example, two unigenes (c127139_g1 and c107331_g1) encoding a multidrug resistance protein (MRP) were not consistent with the pattern of expression in KCd10 vs. KCd0. In addition, for KCd10 vs. BCd10, all unigenes encoding ZIPs, IRT1, HA6t, NRTs, P-type ATPase superfamily (P-ATPase) and copper transporters (COPTs) were down-regulated. The majority of unigenes belonging to ABC transporters (68.4%) and calcium-binding proteins (CMLs, 75.0%) were identified as down-regulated between Baiyewuyueman and Kuishan’aijiaoheiye in Cd treatments. Surprisingly, several critical Cd-related transporters such as natural resistance-associated macrophage proteins (Nramps), zinc induced facilitators (ZIFs) and ABC transporter G family member 35 (PDR8) were similar in the four groups.
Ten transport genes including ZIP2, ZIP3, IRT1, MRP7, MTP3, COPT5, CAX4, HMA2, HMA3 and HMA4 were identified to relevant to the regulation of Cd translocation. The expression of these genes and their biological functions are listed in Table 3. These genes were also used for confirming the reliability of RNA sequencing results by RT-qPCR. As shown in Fig. 6, the relative expression levels of ten genes were significantly higher in BCd10 than in KCd10, and were dramatically up-regulated in BCd10 vs. BCd0. Two strongly up-regulated genes (MTP3 and HMA3) and one significantly down-regulated gene (CAX4) were observed in KCd10 vs. KCd0. Meanwhile, the relative expression of ZIP2, ZIP3, IRT1, COPT5, MRP7, HMA2 and HMA4 were similar between the control and Cd-treatment in Kuishan’aijiaoheiye. The expression patterns from RT-qPCR were similar to those of the RNA-Seq-based gene expression patterns (Fig. 6k). However, three genes including COPT5 (BCd10 vs. BCd0), HMA3 and HMA4 (KCd10 vs. BCd10) did not show consistent expression levels between RT-qPCR and Illumina sequencing data (Fig. 6). The discrepancies may be result from different sensitivity of the two techniques.
Comparative transcriptome analysis of two pak choi cultivars
RNA-Seq approach has been successfully and increasingly used for revealing the mechanisms of Cd resistance and accumulation in plants [21, 12, 26]. According to a recent report, 59,271 unigenes were identified in pak choi, of which 44,539 were functionally annotated and from which many DEGs were obtained . The result provided valuable information for understanding the molecular mechanism of Cd resistance. However, owing to the limitations of genotype and development stages selected, the results could not reflect the overall expression profiles of Cd stress-related genes. In this study, a total of 244,190 unigenes were obtained from eight samples by RNA sequencing and assembly (Table 2), with 142,631 (41.59%) of them successfully annotated by BLAST comparisons with the public databases, which is a greater number than previous reported for pak choi . The results provide a foundation for further investigation of Cd stress mechanisms and identification of novel genes in pak choi. In addition, the large numbers of un-annotated unigenes could be the result of several causes, including novel genes specifically expressed in pak choi roots, the lack of information on the genomes or transcriptomes of pak choi, and other technical or biological biases, such as assembly parameters.
Previous study indicated that the shoot Cd concentration was higher in Baiyewuyueman than in Kuishan’aijiaoheiye under the same Cd condition . In this study, a total of 6827 unigenes were identified as DEGs in four comparisons (Additional file 2: Table S2). Among these DEGs, 368 (368/690, 53.33%) were significantly up-regulated in BCd10 vs. BCd0, and only 99 (99/2733, 3.62%) in KCd10 vs. KCd0. These results suggest that Baiyewuyueman could more effectively activate gene expression under Cd treatment. In addition, 2265 DEGs were specifically regulated by Cd in Kuishan’aijiaoheiye, while only 137 DEGs in Baiyewuyueman (Additional file 2: Table S2). Furthermore, we found that 2239 overlapping DEGs showed similar expression patterns in Kuishan’aijiaoheiye and Baiyewuyueman under Cd stress (Additional file 2: Table S2). These DEGs are not the key genes for regulating Cd tolerance and accumulation in the two cultivars. Besides, 2076 unigenes showed opposing expression patterns in Cd treatment (Additional file 2: Table S2). These results indicate that the two cultivars differed in the molecular mechanisms of Cd response. Additionally, many unigenes were annotated in categories and pathways related to cellular and metabolic processes (Additional file 3: Table S3; Additional file 4: Table S4). These data could be useful for genomic studies and the identification of functional genes in pak choi.
Key genes related to Cd transport in root cells
Cd is mainly transported across the plasma membrane and tonoplast by diverse influx and efflux families of transporters and cation channels in root cells [27, 28]. To date, numerous transporter genes such as ZIP-IRT1 transporters [28, 29], ABC transporters [30, 31], CAXs , OPTs , COPTs  and P1B-ATPase (e.g. HMA2, HMA3 and HMA4) [10, 11, 34], are already known to be involved in Cd transport in plants. Moreover, the role of cation channels (CNGCs and calcium channels) in Cd2+ transport to root cells have been reported [7, 28]. Studies have also demonstrated that multiple transporters play an important role in Cd responses in plants [17, 21, 23]. Xu et al.  identified 24 HM-mediated peptides and transporters as Cd-stressed responsive genes in radish, including members of the CAX2, P-type ATPase and ZIP families. Zhou et al.  identified 63 DEGs belonging to eight GO terms involved in cation transport, which contributed to the genotype differences in Cd accumulating abilities in pak choi. In the current study, many metal transporters and cation channels belonging to various families were also found to be differentially expressed in all pairwise comparisons (Additional file 5: Table S5). Interestingly, 44 unigenes encoding heavy metal transporters, such as ZIP3, ZIP10, IRT1, OPT3, MRP7, PDR9, CAX4, HMA2 and COPT5, showed higher expression in Cd-treated Baiyewuyueman than in Kuishan’aijiaoheiye. However, other heavy metal transporters including HMA (3, 4), PDR (1, 12), MRP2, HKT1 and OPT1, displayed lower expression levels in Baiyewuyueman roots than in Kuishan’aijiaoheiye (Additional file 5: Table S5). In addition, Zhou et al.  found that PDR8 was highly expressed in LAJK compared to HAJS, which may contribute to the low Cd uptake and accumulation in LAJK. OsNRAMP1 and OsNRAMP5 were reported to be Cd transporter involved in the root uptake of Cd from the medium [12,13,14,15]. However, no cultivar difference was detected in the expression of PDR8 and NRAMPs in this study. These results indicate that the mechanisms of Cd2+ transport, involving membrane transporter-related proteins, differed among pak choi cultivars.
Cd transport regulatory networks in pak choi roots
The translocation of Cd in plants is mediated by transporters. Cd generally enters root cells through the transporters of Zn2+, Mn2+, Cu2+ and Fe2+ [7, 35]. Several genes encoding Cd-related transporters have been identified from Arabidopsis [12, 36, 37], rice , wheat  and tobaccos . For example, Cd2+ can enter root cells through ZIP transporters, such as ZIP2, ZIP3 and IRT1 [28, 29]. These transporters have been implicated in mediated Cd influx into root cell across the plasma membrane of root epidermis cells [41, 42]. IRT1 has also been described as an Fe transporter protein that is highly expressed in the cortex cells of Arabidopsis roots . Plants with overexpressed IRT1 accumulated higher levels of Cd and Zn than wild-type plants in Arabidopsis , suggesting that IRT1 also mediates Cd2+ influx into root cells through the plasma membrane of the root epidermis and cortex cells. We found that ZIP2, ZIP3 and IRT1 are expressed at higher levels in Baiyewuyueman than in Kuishan’aijiaoheiye (Fig. 6; Table 3), implying that Baiyewuyueman has higher Cd concentration in roots than Kuishan’aijiaoheiye; however, no significant difference in root Cd concentration was observed between Baiyewuyueman and Kuishan’aijiaoheiye .
The heavy metal P1B-ATPases, such as HMA2 and HMA4 have been identified as efflux transporter genes responsible for the transport of Cd2+ from pericycle parenchyma cells to xylem [14, 15, 44]. These transportors have been localized on the plasma membrane and mainly expressed in the root pericycle cells [39, 45, 46]. In this study, the expression of HMA2 and HMA4 was up-regulated in Baiyewuyueman and down-regulated in Kuishan’aijiaoheiye in Cd treatments compared with the controls (Fig. 6 ; Table 3). Similarly, RT-qPCR analysis revealed that Baiyewuyueman showed higher expression of HMA2 and HMA4 than Kuishan’aijiaoheiye under Cd treatment (Fig. 6). This observation is consistent with previous reports that the expression of HMA4 in Cd-hyperaccumulating species was higher in roots and shoots compared with non-accumulators [47, 48].
Cd is often sequestered in the plant vacuole as Cd-chelate complexes [27, 49, 50]. Several tonoplast-localized transporters, such as CAX4 , HMA3 , MTP3 [52, 53] and MRP7 , are involved in Cd2+ transport into the vacuole. In Arabidopsis, MTP3 has been confirmed to be involved in a Zn translocation, specifically expressed in epidermal and cortex cells of the root , and HMA3 has been identified as involved in Cd translocation, with high expression levels in the guard cells, vascular tissues and root apex . Moreover, the tonoplast-localized AtHMA3 , AtMRP7 , and AtCAX4  have been confirmed to enhance Cd tolerance of plants. RT-qPCR analysis indicates that the four unigenes encoding CAX4, HMA3, MTP3 and MRP7 have higher expression levels in Baiyewuyueman than in Kuishan’aijiaoheiye (Fig. 6a-j), implying the former had higher Cd tolerance than the latter. This result is consistent with our unpublished results that Baiyewuyueman has tolerance to high levels of Cd and rapid growth rate at the seedling stage under Cd treatment. In addition, COPT5 has been identified as a copper export protein located in the tonoplast . In Arabidopsis, Klaumann et al.  showed that COPT5 T-DNA insertion results in decreased copper accumulation in siliques and seeds and increased copper accumulation in the vacuole of root cells. COPT5 has also been confirmed to mediate Cd2+ efflux across the tonoplast of root cells into the cytosolic . Our results show that the expression of COPT5 is higher in Baiyewuyueman than in Kuishan’aijiaoheiye in KCd10 vs. BCd10, and down-regulated in KCd10 vs. KCd0 (Fig. 6 ; Table 3). This finding is consistent with our previous study, implying that shoot Cd concentration was significantly higher in Baiyewuyueman than in Kuishan’aijiaoheiye . Moreover, Cd induced down-regulation for COPT5 and up-regulation for MTP3 and HMA3 in Kuishan’aijiaoheiye, suggesting that these genes might be responsible for enhancing Cd sequestration in the vacuole of root cells in Kuishan’aijiaoheiye. Additionally, many studies have confirmed that the protein systems excluding Cd from the cell cytosol localize to the plasma membrane and tonoplast are mediated by H+-coupled Cd antiport as well as by ATP-energized Cd pumps (e.g. H+-ATPase, V-ATPase and P-ATPase) [7, 28].
Taken together, we proposed a schematic model to depict a regulatory network of root-to-shoot Cd translocation in pak choi (Fig. 7). The low transport rates of Cd into root epidermis cells from soil, and low loading rate of Cd into the xylem sap in the roots, are responsible for the low Cd translocation to shoots in Kuishan’aijiaoheiye. The findings could provide insights into the regulation of Cd translocation in pak choi, and facilitate further manipulation of low Cd accumulation genes in vegetable crops.
The current study identified 6827 unigenes as DEGs in pairwise comparisons. Among them, 690 (BCd10 vs. BCd0), 2733 (KCd10 vs. KCd0), 2919 (KCd0 vs. BCd0) and 3455 (KCd10 vs. BCd10) play significant roles in response to Cd stress. Furthermore, 135 DEGs encoding ion transport related proteins (i.e., ZIPs, P1B-type ATPase and MTPs) were identified. Five plasma membrane-localized transport genes (ZIP2, ZIP3, IRT1, HMA2 and HMA4) and five tonoplast-localized transport genes (CAX4, HMA3, MRP7, MTP3 and COPT5), showing higher expression levels in Baiyewuyueman than in Kuishan’aijiaoheiye, might be responsible for cultivar differences in root-to-shoot Cd translocation in pak choi.
Plant material collection
Pot trials were conducted in a growth chamber with 16 h light at 25 °C and 8 h darkness at 16 °C. Based on the previous work , two pak choi cultivars, Baiyewuyueman (B, high-Cd cultivar) and Kuishan’aijiaoheiye (K, low-Cd cultivar), were used for the pot experiment. In order to avoid damaging the roots at harvest time, sand was used as the culture substrate. The washed sand was spiked with 0 (Cd0) and 10 mg/kg Cd (Cd10) as CdCl2•2.5H2O. The sand (1.2 kg) was filled in each pot (13 cm × 12 cm) and kept moist for 2 weeks prior to the experiment. Seeds were sown directly into pots. The treatments were arranged in a completely randomized design with six replicates (pots).
One week after sowing, the seedlings were thinned to six plants per pot. The plants were fertilized daily with 50 ml of nutrient solution as previously described . Root samples for RNA-seq and RT-qPCR analysis were collected separately at 4 weeks after seedling emergence. Multiple independent biological replicates, each containing a pool of six different plants, were sampled for RNA-Seq (two biological replicates) and RT-qPCR validation (three biological replicates). All samples were immediately frozen in liquid nitrogen and stored at −80 °C.
cDNA library construction and Illumina sequencing
Total RNA of root samples was isolated using Trizol® Reagent (Invitrogen) and purified using the RNeasy Plant Mini kit (Qiagen) according to the manufacturer’s instructions. A total of 3 μg RNA per sample was used for library preparation and sequencing. Eight cDNA libraries named B1_Cd0, B2_Cd0, B1_Cd10, B2_Cd10, K1_Cd0, K2_Cd0, K1_Cd10 and K2_Cd10 were constructed and sequenced as previously described . The library construction and Illumina sequencing were conducted using the Illumina Hiseq™ 2500 platform following the manufacturer’s recommendations at Novogene Bioinformatics Institute (Beijing, China).
The raw reads in fastq format were initially processed with in-house Perl scripts. After filtering the adapter sequences, reads containing ploy-N and low-quality bases, clean reads were obtained. These clean reads were then used for further data analysis.
The left files (read 1 files) from all libraries/samples were pooled into one large left.fq file, and right files (read 2 files) into one large right.fq file. Transcriptome assembly was performed based on the left.fq and right.fq using Trinity  with min_kmer_cov set to 2 by default and all other parameters set as default.
Gene functional annotation and classification
The assembled unigenes were used in BLAST searches against public databases including NCBI non-redundant protein (Nr) and non-redundant nucleotide sequences (Nt), Protein family (Pfam), euKaryotic Ortholog Groups (KOG), Swiss-Prot, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway and Gene Ontology (GO) with an E-value cutoff of 10−5.
Identification and functional enrichment analysis of DEGs
The clean reads were mapped back onto the assembled transcriptome reference sequences by RSEM  (no mismatch) for each sample. Next, the read count for each gene was obtained from the mapping results. Further, the expression level of each gene was normalized to FPKM (fragments per Kilobase per million fragments mapped) based on the number of readcount.
The DEGs were screened as previously described [56, 59]. Prior to differential gene expression analysis, the read counts for each sequenced library were adjusted by edgeR program package through one scaling normalized factor. Next, the average of readcount of the gene from two replicate libraries was calculated as the readcount value of the gene to analyze the differences among four groups: BCd0 (B1_Cd0 and B2_Cd0), BCd10 (B1_Cd10 and B2_Cd10), KCd0 (K1_Cd0 and K2_Cd0) and KCd10 (K1_Cd10 and K2_Cd10). The DESeq R package (1.10.1) was used to identify DEGs between two groups according to the method described by Anders & Huber . DESeq provides statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution. The resulting P-values were adjusted using the Benjamini and Hochberg’s approach for controlling the false discovery rate. A unigene found by DESeq with an adjusted P-value <0.05 was considered significantly differentially expressed.
Functional enrichment analysis was performed by the GOseq R packages based Wallenius non-central hyper-geometric distribution (corrected P-value <0.05)  and KOBAS (2.0) software (corrected P-value <0.05) . Cluster analysis of gene expression patterns was performed with cluster 3.0 and Java Treeview software.
Total RNA was isolated as described above. First strand cDNA fragments were synthesized using the PrimeScript® RT reagent kit (Takara, Dalian, China). RT-qPCR was performed on an ABI 7300 (Applied Biosystems, Foster City, CA, USA) using a SYBR Premix EX Taq kit (Takara) in a 20 μl reaction mixtures, containing 2 μl of diluted cDNA, 0.2 μM forward and reverse primer, and 10 μl 2 × SYBR Green PCR Master Mix. The PCR reaction protocol was 95 °C for 30 s, 40 cycles of 95 °C for 5 s and 60 °C for 30 s. The fluorescence was measured via a 65–95 °C melting curve. All reactions were performed in triplicate with each replicate measured three times; the relative expression level of the selected unigenes using the Actin gene as the internal control gene was calculated using ratio = 2−ΔΔCτ. The specific primers for RT-qPCR were designed using Beacon Designer 7.0 software (Premier Biosoft International, USA) and are listed in Additional file 6: Table S6.
Cyclic nucleotide gated channel
Differentially expressed gene
Cadmium/zinc/copper-transporting ATPase (heavy metal-associated domain)
Kyoto encyclopedia of genes and genomes
Multidrug resistance protein
Metal tolerance protein
Natural resistance-associated macrophage protein
Pleiotropic drug resistance-type ABC transporter protein
Reverse-transcription quantitative PCR
Zinc induced facilitator
Zinc-regulated transporter/Iron-regulated transporter-like Protein
Wang M, Zou J, Duan X, Jiang W, Liu D. Cadmium accumulation and its effects on metal uptake in maize (Zea mays L.). Bioresour Technol. 2007;98(1):82–8.
Xu L, Wang Y, Zhai L, Xu Y, Wang L, Zhu X, et al. Genome-wide identification and characterization of cadmium-responsive microRNAs and their target genes in radish (Raphanus sativus L.) roots. J Exp Bot. 2013;64(14):4271–87.
Chen X, Wang J, Shi Y, Zhao MQ, Chi GY. Effects of cadmium on growth and photosynthetic activities in pakchoi and mustard. Bot Stud. 2011;52:41–6.
Clemens S, Aarts MG, Thomine S, Verbruggen N. Plant science: the key to preventing slow cadmium poisoning. Trends Plant Sci. 2013;18(2):92–9.
Sherameti I, Varma A. Heavy metal contamination of soils: monitoring and remediation. 4th ed. New York: Springer; 2015.
Xiong J, An L, Lu H, Zhu C. Exogenous nitric oxide enhances cadmium tolerance of rice by increasing pectin and hemicellulose contents in root cell wall. Planta. 2009;230(4):755–65.
Lux A, Martinka M, Vaculík M, White PJ. Root responses to cadmium in the rhizosphere: a review. J Exp Bot. 2011;62(1):21–37.
Li T, Tao Q, Shohag MJI, Yang X, Sparks DL, Liang Y. Root cell wall polysaccharides are involved in cadmium hyperaccumulation in Sedum alfredii. Plant Soil. 2015;389(1–2):387–99.
Satoh-Nagasawa N, Mori M, Nakazawa N, Kawamoto T, Nagato Y, Sakurai K, et al. Mutations in rice (Oryza sativa) heavy metal ATPase 2 (OsHMA2) restrict the translocation of zinc and cadmium. Plant Cell Physiol. 2012;53(1):213–24.
Miyadate H, Adachi S, Hiraizumi A, Tezuka K, Nakazawa N, Kawamoto T, et al. OsHMA3, a P1B-type of ATPase affects root-to-shoot cadmium translocation in rice by mediating efflux into vacuoles. New Phytol. 2011;189:190–9.
Takahashi R, Bashir K, Ishimaru Y, Nishizawa NK, Nakanishi H. The role of heavy-metal ATPases, HMAs, in zinc and cadmium transport in rice. Plant Signal Behav. 2012;7(12):1605–7.
Takahashi R, Ishimaru Y, Senoura T, Shimo H, Ishikawa S, Arao T, et al. The OsNramp1 iron transporter is involved in Cd accumulation in rice. J Exp Bot. 2011;62(14):4843–50.
Sasaki A, Yamaji N, Yokosho K, Ma JF. Nramp5 is a major transporter responsible for manganese and cadmium uptake in rice. Plant Cell. 2012;24(5):2155–67.
Ishimaru Y, Takahashi R, Bashir K, Shimo H, Senoura T, Sugimoto K, et al. Characterizing the role of rice NRAMP5 in manganese, iron and cadmium transport. Sci Rep. 2012;2(6071):286.
Ishikawa S, Ishimaru Y, Igura M, Kuramata M, Abe T, Senoura T, et al.Ion-beam irradiation, gene identification, and marker-assisted breeding in the development of low-cadmium rice. Proc Natl Acad Sci U S A. 2012;109(47):19166–71.
Ward JA, Ponnala L, Weber CA. Strategies for transcriptome analysis in nonmodel plants. Am J Bot. 2012;99(2):267–76.
Xu L, Wang Y, Liu W, Wang J, Zhu X, Zhang K, et al. De novo sequencing of root transcriptome reveals complex cadmium-responsive regulatory networks in radish (Raphanus sativus L.). Plant Sci. 2015;236:313–23.
Herbette S, Taconnat L, Hugouvieux V, Piette L, Magniette ML, Cuine S, et al. Genome-wide transcriptome profiling of the early cadmium response of Arabidopsis roots and shoots. Biochimie. 2006;88(11):1751–65.
Liu T, Zhu S, Tang Q, Tang S. Genome-wide transcriptomic profiling of ramie (Boehmeria nivea L. Gaud) in response to cadmium stress. Gene. 2015;558(1):131–7.
Gao J, Luo M, Zhu Y, He Y, Wang Q, Zhang C. Transcriptome sequencing and differential gene expression analysis in Viola yedoensis Makino (Fam. Violaceae) responsive to cadmium (Cd) pollution. Biochem Biophys Res Commun. 2015;459(1):60–5.
Xu J, Sun JH, Du LG, Liu XJ. Comparative transcriptome analysis of cadmium responses in Solanum nigrum and Solanum torvum. New Phytol. 2012;196(1):110–24.
Cao FB, Chen F, Sun H, Zhang G, Chen ZH, Wu F. Genome-wide transcriptome and functional analysis of two contrasting genotypes reveals key genes for cadmium tolerance in barley. BMC Genomics. 2014;15(1):611.
Zhou Q, Guo JJ, He CT, Shen C, Huang YY, Chen JX, et al. Comparative transcriptome analysis between low- and high-cadmium-accumulating genotypes of pakchoi (Brassica chinensis L.) in response to cadmium stress. Environ Sci Technol. 2016;50(12):6485–94.
Chen Y, Li TQ, Han X, Ding ZL, Yang XE, Jin YF. Cadmium accumulation in different pakchoi cultivars and screening for pollution-safe cultivars. J Zhejiang Univ Sci B. 2012;13(6):494–502.
Xia S, Deng R, Zhang Z, Liu C, Shi G. Variations in the accumulation and translocation of cadmium among pak choi cultivars as related to root morphology. Environ Sci Pollut Res. 2016;23(10):9832–42.
Milner MJ, Mitani-Ueno N, Yamaji N, Yokosho K, Craft E, Fei Z, et al. Root and shoot transcriptome analysis of two ecotypes of Noccaea caerulescens uncovers the role of NcNramp1 in Cd hyperaccumulation. Plant J. 2014;78(3):398–410.
Gallego SM, Pena LB, Barcia RA, Azpilicueta CE, Iannone MF, Rosales EP, et al. Unravelling cadmium toxicity and tolerance in plants: insight into regulatory mechanisms. Environ Exp Bot. 2012;83:33–46.
Asgher M, Khan MIR, Anjum NA, Khan NA. Minimising toxicity of cadmium in plants-role of plant growth regulators. Protoplasma. 2015;252(2):399–413.
Eide D, Broderius M, Fett J, Guerinot ML. A novel iron-regulated metal transporter from plants identified by functional expression in yeast. Proc Natl Acad Sci. 1996;93(11):5624–8.
Park J, Song WY, Ko D, Eom Y, Hansen TH, Schiller M, et al. The phytochelatin transporters AtABCC1 and AtABCC2 mediate tolerance to cadmium and mercury. Plant J. 2012;69(2):278–88.
Wojas S, Hennig J, Plaza S, Geisler M, Siemianowski O, Skłodowska A, et al. Ectopic expression of Arabidopsis ABC transporter MRP7 modifies cadmium root-to-shoot transport and accumulation. Environ Pollut. 2009;157(10):2781–9.
Wu Q, Shigaki T, Williams KA, Han JS, Kim CK, Hirschi KD, et al. Expression of an Arabidopsis Ca2+/H+ antiporter CAX1 variant in petunia enhances cadmium tolerance and accumulation. J Plant Physiol. 2011;168(2):167–73.
Klaumann S, Nickolaus SD, Fürst SH, Starck S, Schneider S, Ekkehard Neuhaus H, et al. The tonoplast copper transporter COPT5 acts as an exporter and is required for interorgan allocation of copper in Arabidopsis thaliana. New Phytol. 2011;192(2):393–404.
Verret F, Gravot A, Auroy P, Leonhardt N, David P, Nussaume L, et al. Overexpression of AtHMA4 enhances root-to-shoot translocation of zinc and cadmium and plant metal tolerance. FEBS Lett. 2004;576(3):306–12.
Komal T, Mustafa M, Ali Z, Kazi AG. Heavy metal uptake and transport in plants. In: Sherameti I, Varma A, editors. Heavy metal contamination of soils. New York: Springer International Publishing; 2015. p. 181–94.
Carrió-Seguí A, Garcia-Molina A, Sanz A, Peñarrubia L. Defective copper transport in the copt5 mutant affects cadmium tolerance. Plant Cell Physiol. 2015;56(3):442–54.
Morel M, Crouzet J, Gravot A, Auroy P, Leonhardt N, Vavasseur A, et al. AtHMA3, a P1B-ATPase allowing Cd/Zn/Co/Pb vacuolar storage in Arabidopsis. Plant Physiol. 2009;149(2):894–904.
Takahashi R, Ishimaru Y, Shimo H, Ogo Y, Senoura T, Nishizawa NK, et al. The OsHMA2 transporter is involved in root-to-shoot translocation of Zn and Cd in rice. Plant Cell Environ. 2012;35(11):1948–57.
Tan J, Wang J, Chai T, Zhang Y, Feng S, Li Y, et al. Functional analyses of TaHMA2, a P(1B)-type ATPase in wheat. Plant Biotechnol J. 2013;11(4):420–31.
Koren’kov V, Park S, Cheng NH, Sreevidya C, Lachmansingh J, Morris J, et al. Enhanced Cd2+-selective root-tonoplast-transport in tobaccos expressing Arabidopsis cation exchangers. Planta. 2007;225(2):403–11.
Milner MJ, Seamon J, Craft E, Kochian LV. Transport properties of members of the ZIP family in plants and their role in Zn and Mn homeostasis. J Exp Bot. 2013;64(1):369–81.
Vert G, Grotz N, Dédaldéchamp F, Gaymard F, Guerinot ML, Briat J-F, et al. IRT1, an Arabidopsis transporter essential for iron uptake from the soil and for plant growth. Plant Cell. 2002;14(6):1223–33.
Connolly EL, Fett JP, Guerinot ML. Expression of the IRT1 metal transporter is controlled by metals at the levels of transcript and protein accumulation. Plant Cell. 2002;14(6):1347–57.
Mills RF, Francini A, Ferreira da Rocha PS, Baccarini PJ, Aylett M, Krijger GC, et al. The plant P1B-type ATPase AtHMA4 transports Zn and Cd and plays a role in detoxification of transition metals supplied at elevated levels. FEBS Lett. 2005;579(3):783–91.
Park W, Ahn SJ. How do heavy metal ATPases contribute to hyperaccumulation?. J Plant Nutr Soil Sci. 2014;177(2):121–7.
Ueno D, Milner MJ, Yamaji N, Yokosho K, Koyama E, Zambrano MC, et al. Elevated expression of TcHMA3 plays a key role in the extreme Cd tolerance in a Cd-hyperaccumulating ecotype of Thlaspi caerulescens. Plant J. 2011;66(5):852–62.
Bernard C, Roosens N, Czernic P, Lebrun M, Verbruggen N. A novel CPx-ATPase from the cadmium hyperaccumulator Thlaspi caerulescens. FEBS Lett. 2004;569(1–3):140–8.
Papoyan A, Kochian LV. Identification of Thlaspi caerulescens genes that may be involved in heavy metal hyperaccumulation and tolerance. Characterization of a novel heavy metal transporting ATPase. Plant Physiol. 2004;136(3):3814–23.
Liu W, Zhang X, Liang L, Chen C, Wei S, Zhou Q. Phytochelatin and oxidative stress under heavy metal stress tolerance in plants. In: Gupta DK, Palma JM, Corpas FJ, editors. Reactive oxygen species and oxidative damage in plants under stress. New York: Springer; 2015. p. 191–217.
Shi YZ, Zhu XF, Wan JX, Li GX, Zheng SJ. Glucose alleviates cadmium toxicity by increasing cadmium fixation in root cell wall and sequestration into vacuole in Arabidopsis. J Integr Plant Biol. 2015;57(10):830–7.
Korenkov V, King B, Hirschi K, Wagner GJ. Root-selective expression of AtCAX4 and AtCAX2 results in reduced lamina cadmium in field-grown Nicotiana tabacum L. Plant Biotechnol J. 2009;7(3):219–26.
Sharma SS, Dietz KJ, Mimura T. Vacuolar compartmentalization as indispensable component of heavy metal detoxification in plants. Plant Cell Environ. 2016;39(5):1112–26.
Arrivault S, Senger T, Krämer U. The Arabidopsis metal tolerance protein AtMTP3 maintains metal homeostasis by mediating Zn exclusion from the shoot under Fe deficiency and Zn oversupply. Plant J. 2006;46(5):861–79.
Gravot A, Lieutaud A, Verret F, Auroy P, Vavasseur A, Richaud P. AtHMA3, a plant P1B-ATPase, functions as a Cd/Pb transporter in yeast. FEBS Lett. 2004;561(1–3):22–8.
Garcia-Molina A, Andrés-Colás N, Perea-García A, del Valle-Tascón S, Peñarrubia L, Puig S. The intracellular Arabidopsis COPT5 transport protein is required for photosynthetic electron transport under severe copper deficiency. Plant J. 2011;65(6):848–60.
Hu L, Li H, Chen L, Lou Y, Amombo E, Fu J. RNA-seq for gene identification and transcript profiling in relation to root growth of bermudagrass (Cynodon dactylon) under salinity stress. BMC Genomics. 2015;16(1):575.
Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29(7):644–52.
Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12(1):93–9.
Liu MM, Xing YM, Zhang DW, Guo SX. Transcriptome analysis of genes involved in defence response in Polyporus umbellatus with Armillaria mellea infection. Sci Rep. 2015;5:16075.
Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11(10):R106.
Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11(2):R14.
Mao X, Cai T, Olyarchuk JG, Wei L. Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics. 2005;21(19):3787–93.
This work was in part supported by grants from the Natural Science Foundation of China (Nos. 31370515 and 31671599), Natural Science Foundation of Anhui Province of China (1708085MC55) and the Science Research Foundation of the Education Bureau of Anhui Province of China (KJ2016B014).
Availability of data and materials
The RNA sequence dataset supporting the results of this article is available in the repository of NCBI Sequence Read Archive (SRA) with the GenBank accession No.: SRX2442157, SRX2442158, SRX2442159 and SRX2442163.
Ethics approval and consent to participate
Consent for publication
All the authors have signed the consent form.
The authors declare that they have no competing interests for this research.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
KEGG pathways of the assembled unigenes in pak choi. (XLSX 239 kb)
Complete lists of differentially expressed genes in four groups based on pairwise comparisons. (XLSX 6819 kb)
The enriched GO terms of up and down-regulated DEGs in four comparisons. (XLSX 153 kb)
The significantly enriched pathways of up and down-regulated DEGs in four groups based on pairwise comparison. (XLSX 16 kb)
The DEGs involved in transport activities in two pak choi cultivars under Cd stress. (XLSX 84 kb)
The specific primer sequences of ten DEGs validated by RT-qPCR analysis. (XLSX 11 kb)
About this article
Cite this article
Yu, R., Li, D., Du, X. et al. Comparative transcriptome analysis reveals key cadmium transport-related genes in roots of two pak choi (Brassica rapa L. ssp. chinensis) cultivars. BMC Genomics 18, 587 (2017). https://doi.org/10.1186/s12864-017-3973-2
- Pak choi
- Cd stress