Comparative transcriptomics provides novel insights into the mechanisms of selenium tolerance in the hyperaccumulator plant Cardamine hupingshanensis

Selenium (Se) is an essential mineral element for animals and humans. Cardamine hupingshanensis (Brassicaceae), found in the Wuling mountain area of China, has been identified as a novel Se hyperaccumulator plant. However, the mechanism for selenium tolerance in Cardamine plants remains unknown. In this study, two cDNA libraries were constructed from seedlings of C. hupingshanensis treated with selenite. Approximately 100 million clean sequencing reads were de novo assembled into 48,989 unigenes, of which 39,579 and 33,510 were expressed in the roots and leaves, respectively. Biological pathways and candidate genes involved in selenium tolerance mechanisms were identified. Differential expression analysis identified 25 genes located in four pathways that were significantly responsive to selenite in C. hupingshanensis seedlings. The results of RNA sequencing (RNA-Seq) and quantitative real-time PCR (RT-qPCR) confirmed that storage function, oxidation, transamination and selenation play very important roles in the selenium tolerance in C. hupingshanensis. Furthermore, a different degradation pathway synthesizing malformed or deformed selenoproteins increased selenium tolerance at different selenite concentrations. This study provides novel insights into the mechanisms of selenium tolerance in a hyperaccumulator plant, and should serve as a rich gene resource for C. hupingshanensis.

Selenium (Se) is an essential trace element for animals and humans that can be acquired from plant accumulators growing in seleniferous soil. According to tolerance and accumulation quantities of Se, plants can be categorized into three groups: <100 mg Se kg −1 , 100-1000 mg Se kg −1 and 1000-15000 mg Se kg −1 . Plants which can tolerate or accumulate Se at a concentration of 1000-15000 mg Se kg −1 are called Se hyperaccumulators 1 . Most species known to hyperaccumulate Se belong to the Fabaceae family. The ability of hyperaccumulation of Se in plants has evolved several times within the Asteraceae, Brassicaceae and Fabaceae 1 . Astragalus bisulcatus, Stanleya pinnata and Symphyotrichum ericoides are the most widely studied Se hyperaccumulators [2][3][4][5][6][7] .
Cardamine hupingshanensis is a novel hyperaccumulator plant found in the Wuling mountain area. Bai et al. 8 found that C. hupingshanensis is primarily distributed in Hunan province in China, at 800-1400 m. However, we found it grows where there is a cloudy slope or valley with coal gangue and running water in the city of Enshi as well as in the counties of Xuan' en, Changyang and Wufeng in Hubei province at 800-1900 m. Yuan et al. (2013) and Shao et al. (2014) measured concentrations of total Se by hydride generation-atomic fluorescence spectrometry (HG-AFS) and HPLC-ICP-MS, respectively. These studies showed that C. hupingshanensis could accumulate Se in excess of 1400 mg Se kg −1 of dry matter in all tissues of seedlings, with most not exceeding 4000 mg Se kg −1 of dry matter in roots 9,10 .

Results
Transcriptome characteristics in C. hupingshanensis. RNA samples from leaves and roots of C. hupingshanensis were prepared for library construction and subsequently sequenced on the Illumina HiSeq. 2500 platform. We obtained a total of 54,765,658 and 50,352,860 raw paired-end reads in leaves and roots, respectively. All sequencing data were deposited in the NCBI database and can be accessed with the Sequence Read Archive (SRA) number of SRP097726. After quality analysis and data filtering, 52,019,342 and 48,078,676 clean reads were retained with Q20 values of 99.0% and 98.9% and GC contents of 45.3% and 47.3% in the leaves and roots, respectively. We performed a de novo transcript assembly using these paired-end data to obtain transcript sequences. A total of 78,471 transcripts (average GC content of 41.18%), including 48,989 unigenes, were assembled with a total length of 86,620,844 bp (average transcript length of 1,104 bp). The size and copy distribution of the transcripts are displayed in Fig. 2a and b. The transcript abundance analyzed by bowtie (version 2.23) and RSEM (version 1.2.15) showed that 39,579 and 33,510 transcripts were expressed in the roots and leaves, respectively (displayed in Fig. 2c and d). The function of each unigene set in C. hupingshanensis was then annotated by Trinotate (version r20131110) based on homologies to putative or known sequences available in public databases (Table 1). In addition, a gene ontology (GO) analysis, which is a major bioinformatic approach utilizing to represent properties of gene and gene products across all species, was then carried out on the putative proteins. All unigenes were annotated using three ontologies, including biological process (BP), molecular function (MF) and cell component (CC) (Fig. 3). There was a total of 51 different sublevels narrowed down to form the three ontologies. According to the explanations in the non-redundant protein (NR) and the Pfam databases, 48,989 unigenes properly fit into one or more ontologies.
There were 23 BP subcategories with 40,081 unigenes, 13 CC subcategories with 36,783 unigenes and 15 MF subcategories with 20,031 unigenes. Consistent with findings in other plants, the metabolic process, cellular process and single cell process ontologies in BP were the top three gene ontology terms, with 9,167, 8,555 and 5,581 unigenes, respectively 22,23 . Cell, cell part and organelle terms in CC were the top three classes with 10,818, 10,818 and 7,638 unigenes, respectively. Catalytic activity, nucleic acid binding and transcription factor activity in MF were the top three GO terms with 9,009, 7,543 and 1,313 unigenes, respectively. The COG function classification of the C. hupingshanensis unigenes is displayed in Fig. 4. Overall, 14,417 of 48,989 unigenes matched to the COG database were clustered into 24 functional clusters. According to the number of genes, the most significant cluster was the general function prediction only cluster (2,463, 17.08%), followed by the nucleotide transporter and  To understand the functions and products of unigenes in putative metabolic pathways, the Kyoto Encyclopedia of Genes and Genomes (KEGG) was used to systematically analyze all unigenes in C. hupingshanensis (Fig. 5). A total of 5,196 unigenes obtained in this study were classified into five branches: cellular processes, environmental information processing, genetic information processing, metabolism, and organismal systems. More than half of those aligned with KEGG transcripts were classified into metabolism (59.37%), and 22.19% were classified into environmental information processing. The other highly represented pathways included the global and overview map (1,180,22 General DEGs at all selenium concentrations. To understand the mechanisms of selenium tolerance in the hyperaccumulator plant C. hupingshanensis, another eighteen libraries (one sample including 3 biological replicates, nine libraries for roots and nine libraries for leaves) of seedlings were constructed to identify differentially expressed genes (DEGs) between the control and low Se treatment (100 μg Se/L, slightly higher than the concentration of Se in the water of the high-Se area, treated for 24 hours) and between the control and high Se  treatment (80,000 μg Se/L, a stress concentration, treated 24 hours). Using a probability value of more than 0.8 and a minimal FPKM value of 3 (FPKM values used to define up-and down-regulated genes between treatment and control), the overlapping parts and exclusive sections of 670 unigenes from the four comparative groups are shown in Fig. 6. There were 50, 181, 95, and 128 unigenes annotated exclusively in these four groups. There were 4 unigenes that aggregated into a collection from the four groups (Supplementary Table S1). There were 43 and 50 unigenes transcribed in both the roots and leaves at two selenium concentrations, respectively. The roots and leaves shared 9 and 13 unigenes at low and high selenium concentrations, respectively.
Using a probability value of more than 0.8 and a minimal FPKM value of 10 (FPKM values used to define upand down-regulated genes between treatment and control), 31 and 64 genes were associated with Se response in the root tissue, and 30 (23 genes up-regulated and 7 genes down-regulated) and 103 (42 genes up and 61 genes down) genes were associated with Se response in leaf tissue of seedlings treated with low and high concentrations  of Se, respectively. At both concentrations of Se treatment, 6 annotated genes were up-regulated and 1 was down-regulated in root tissue, and 8 genes (including 3 un-annotated genes) were up-regulated and 6 (including 1 un-annotated gene) were down-regulated in the leaf tissue (Table 2).
In this study, the expression of sulfite oxidase (SOX) gene in the root was up-regulated when selenite was added to the culture solution of the C. hupingshanensis seedlings ( Table 2 and Fig. 7A), suggesting that selenite may be converted to selenate first and then the selenate continued through metabolism.
Vacuole is of great importance because of its storage function, which was demonstrated by the up-regulation of genes in metabolic and transport pathways. Based on the data from RNA-Seq and RT-qPCR (Table 2, Fig. 7A and B), the changes in genes of the glutathione S-transferase family and the C subfamily of the ATP-binding cassette transporters (ABCC) provided direct evidence. The glutathione S-transferase family gene GST-u4 in leaves, which promotes glutathione-chelated selenate to form phytochelatins (PCs) 24 , was up-regulated. The expression of ABCC2 was significantly up-regulated in leaves of C. hupingshanensis seedlings when the seedlings were treated with selenite. At the same time, there were four up-regulated genes associated with metal ion binding.
The results of RNA-Seq and RT-qPCR (Table 3, Fig. 7C) showed that the oxidation and transamination of SeCys might be two important processes for Se detoxification in the roots, but the results of RT-qPCR (Fig. 7D) suggested that conversion of SeCys to SeMet was also an important process for Se detoxification in the leaves. The transamination of SeCys or its oxides, such as L-cysteate and L-cysteine-sulfinate, has not been studied extensively in plants. The aspartate aminotransferase (chloroplastic, Aat), a pyridoxal phosphate dependent amino acid aminotransferase, which also located in cysteine and methionine metabolism pathway, could catalyze SeCys, L-cysteate and L-cystiene-sulfinate deaminize [25][26][27] . It was found that the expression of Aat and pyridoxine 4-dehydrogenase (PLR1) were up-regulated in the roots. Cystathionine gamma-synthase (CγS), cystathionine beta-lyase (CBL) and 5-methyltetrahydropteroyl-triglutamate-homocysteine methyltransferase (MET), which are three key enzymes in the methionine biosynthesis pathway, were found to be up-regulated by RT-qPCR at all concentrations of selenium treatment in leaves (Fig. 7D).
Selenation which is similar to sulfation in plants, is another important pathway for detoxification of selenium because it can reduce SeCys biosynthesis. The aryl sulfotransferase enzyme (SULT1A, EC 2.8.2.1) is a member of the sulfotransferase (SOT, EC 2.8.2.-) protein family. This enzyme can transfer a sulfate group from the donor 3′-phosphoadenosine 5′-phosphosulfate (PAPS) to phenolic sulfate esters to a phenolic acceptor substrate in the sulfur and glucosinolate metabolism pathway. According to the KEGG analysis, this enzyme plays a key role in the selenation pathway 28 . However, aryl sulfotransferase is present in animals and prokaryotes, while aryl sulfotransferase is found only in R. communis L. [28][29][30] . The expression of SULT1A increased 4.94 (100 μg Se/L)-and 7.17 (80,000 μg Se/L)-fold in the leaves ( Table 2 and Fig. 7F). That is, the response of selenation is a common pathway in leaves of C. hupingshanensis when seedlings are treated with selenium. Therefore, we can conclude that selenation is the more common approach of tolerating selenium stress in leaves of C. hupingshanensis seedlings.  The other genes that significantly responded to Se are involved in multiple physiological processes. The gene for transcription factor LBD 16 of the LBD family, which plays crucial roles in diverse growth and development processes, including the establishment and maintenance of the developmental boundary of lateral organs 31 , was up-regulated in the roots, indicating that selenium could promote lateral root formation in plants. AtGoLS3 is a member of the galactinol synthase (GolS) family, which initiates the biosynthesis of raffinose oligosaccharides (RFO), may act as an osmoprotectant in drought stress tolerance through UDP-galactose 32,33 . According to results of RNA-Seq, a novel possibility is that AtGoLS3 is induced not only by cold stress but also by selenium stress. The cellular-localized cold-inducible RNA-binding protein, also called the glycine-rich RNA-binding protein, which appears to be involved in the adaptation of abiotic and biotic stress 34 , was another significantly up-regulated gene that strongly responded to selenium in the leaves. Another significantly up-regulated gene that responded to selenium was that coding for the dormancy/auxin-associated protein, which is involved in growth suppression in bud and hypocotyl tissues and whose expression increases in response to abiotic or biotic factors 35 . However, xyloglucan xyloglucosyl transferase, also called xyloglucan endo-transglycosylase (XET), which is involved in wall-loosening, wall-strengthening process, gravitropic responses and the incorporation of nascent xyloglucan into the wall during biosynthesis 36 , was significantly down-regulated at all concentrations of selenium. The down-regulated regulatory beta subunit was one of three subunits of 5′-AMP-activated protein kinase (AMPK), which plays a role in maintaining cellular energy homeostasis 37 . Selenium also affected the process of phagolysosome formation by down-regulation of tubulin beta units of the microtubule. Special DEGs at low concentrations of Se treatment. The same criteria described above were used to analyze special DEGs in the seedlings of C. hupingshanensis treated with low Se concentrations (100 μg Se/L). There were 9 up-regulated genes and 5 down-regulated genes that were found to be specifically expressed in the roots, and 7 up-regulated genes and 2 down-regulated genes were found to be specifically expressed in the leaves. The low Se concentration treatment not only significantly increased the storage function of the vacuole as well as selenation and transamination but also invoked protein degradation and other physiological responses (Table 3 and Fig. 7G). E3 ubiquitin-protein ligase RNF13 and MUL1, which are two important members of the ubiquitin-proteasome pathway (UPP), were significantly up-regulated in the roots under low Se treatment. RING finger protein 13 (RNF13) is a member of the largest family of ubiquitin ligases in eukaryotes and is an ER/Golgi membrane-associated E3 ubiquitin ligase that has been identified as a novel RING-based ubiquitin ligase 38,39 . Mitochondrial E3 ubiquitin protein ligase 1 (MUL1), which is localized to the mitochondria, is a crucial moderator of retinoic acid-inducible-gene I (RIG-I) signaling 40 . Cysteine protease RD19A, which is an important enzyme of proteolysis involved in cellular protein catabolic processes and responses to osmotic and salt stress that is located mainly in the vacuole and lysosome 41 , was also up-regulated (Table 3 and Fig. 7G). The low concentration of Se treatment significantly affected carbohydrate metabolism by up-regulating phosphoenolpyruvate carboxykinase and the small chain of ribulose-bisphosphate carboxylase in the roots and beta-amylase and beta-glucosidase in the leaves. Some genes related to the circadian clock, such as GIGANTEA (GI) 42,43 and pseudo-response regulator 5 (PRR5) 44 as well as those of the two-component response regulator ARR-B family (PCL1) 45,46 , were up-regulated, which may suggest effects of selenium on the circadian clock in leaves of C. hupingshanensis. The two least up-regulated genes were saposin and ATP-dependent RNA helicase DDX5/DBP2, and DDX5/DBP2 act in the process of nonsense-mediated mRNA decay and ribosome biogenesis through rRNA. A few genes including brassinosteroid insensitive 1-associated receptor kinase 1, time for coffee, endonuclease 2 and ribulose-bisphosphate carboxylase small chain were inhibited by low Se treatment in the roots. Only golgin subfamily A member 6-like protein 22 and a putative transcription factor were down-regulated in the leaves.
Special DEGs at high concentrations of Se treatment. The same criteria described above were used to analyze special DEGs in C. hupingshanensis seedlings treated by a high concentration of Se (80,000 μg Se/L). Compared with the low Se concentration treatment, there were more genes with expression changes of various physiological processes and cellular functions: 23 genes were up-regulated and 31 genes were down-regulated specifically in the roots, and 22 genes were up-regulated and 12 genes were down-regulated specifically in the leaves ( Table 4).
The plant C. hupingshanensis, a novel selenium hyperaccumulator, had distinctive reactions to high concentrations of selenium (Table 4 and Fig. 7D). The first reaction was repression of selenium uptake through down-regulating the expression of sulfate transporter 1.2 (Sultr1;2) which was a key protein involved in sulfate and selenate transport and expressed mainly in the root cortex, the root tip and lateral roots 47,48 . The second change was repressed reduction of selenate in roots by down-regulating the expression of adenylyl-sulfate reductase (glutathione, APR1), a critical enzyme catalyzing reduction of adenosine 5′ phoshposulfate (APS) or phoshposelenate (APSe), in which C and N terminal domains had a GRX and TRX like function respectively 49 . The last distinctive change comes from the up-regulated expression of aryl sulfotransferase (SULT1A), which is closely connected with sulfation or selenation. These results indicated that the flux of selenate on APS was converted to  PAPSe and used for selenation, but was not reduced to selenide or combined into SeCys and selenoprotien in the root when C. hupingshanensis was treated with high concentrations of selenium 50 . The visible effects of redox homeostasis when treated with high concentrations of selenium were through the regulation of thioredoxin (TRX) and glutaredoxin (GRX) which were involved in detoxifying reactive oxygen species (ROS) during stress responses and determination protein thiol/disulfide status, and played key roles in the maintenance of cellular redox homeostasis through the sensing and reducing equivalents to a large number of target proteins, such as reductases, peroxidases, transcription factors, metabolic enzymes of glycolysis, and photosynthesis or through structural modifications of target proteins 51,52 . The genes of thioredoxin 1 (Trx 1) and glutaredoxin C-10 (GrxC10) from the roots were up-regulated (Table 4 and Fig. 7D) when the seedlings of C. hupingshanensis were treated by high concentrations of selenite. Simultaneously, peroxidase 43 and 67 which were the important target proteins of TRX, and GRX were up-regulated in the roots 53 .
Regarding photosynthesis, high selenium concentration suppressed the expression of genes involved in light and dark reaction. CP43, one of the components of the core complex of photosystem II (PSII) which binds chlorophyll and helps catalyze the primary light-induced photochemical processes of PSII 54 , was down-regulated when treated with high concentrations of selenium in the leaves. PsaA and PsaB, which bind P700 and are the primary electron donor of photosystem I (PSI), as well as the electron acceptors A0, A1 and FX 55 , were also  down-regulated when treated with high concentrations of selenium in the leaves. So, not only the light harvesting process located in PSII but also the transferring of electron process located in PSI of the light reaction of photosynthesis were suppressed when treated with high concentrations of selenium. On the other hand, Rubisco is the key enzyme complex in dark reaction and catalyzes two reactions: the carboxylation of D-ribulose 1,5-bisphosphate, the primary event in carbon dioxide fixation, as well as the oxidative fragmentation of the pentose substrate, but the L subunit of Rubisco was down-regulated when treated with high concentrations of selenium in the leaves. The carbon fixation in dark reaction was also inhibited by high concentrations of selenium through down-regulating the gene expression of the L subunit of Rubisco. In addition, phytochrome-interacting factor 3 (PIF3) [56][57][58] , which is a basic helix-loop-helix (bHLH) transcription factor closely related to the switch between skotomorphogenesis and photomorphogenesis, was up-regulated in the roots when the seedlings were exposed to light, but the function of PIF3 decided by its state phosphorylation, is still unclear.
The root growth and development of C. hupingshanensis seedlings were affected predominantly by high concentrations of selenium. The first change comes from the genes involved in lignin biosynthesis. Five members of class III perooxidases were changed, which play critical roles in lignin biosynthesis, reduction of hydrogen peroxide, auxin and secondary metabolism 53 . The genes of peroxidase 43 and 67 were up-regulated, and peroxidase 3, 11a and 56 were down-regulated in the roots. At the same time, the genes encoding ferulate-5-hydroxylase and coniferyl-aldehyde dehydrogenase in the phenylpropanoid biosynthesis pathway associated with the production of precursors for lignin biosynthesis were down-regulated 59,60 . Additionally, the down-regulation of respiratory burst oxidase homologs (rbohs) decreases the production of superoxide 61 in the roots. The second change comes from the genes involved in the process of the development of roots. Gibberellin 2-oxidase (GA2ox), which plays very important roles in plant growth and development and can alter expression of lignin biosynthesis-related genes to reduce biomass accumulation and lignification 62 , was up-regulated; it is also involved in resistance to high-salinity stress 63 . All the above changes may indicate that selenium can affect the rigidity and strength of roots. The up-regulation of transcription factor LBD 16, auxin-responsive protein IAA (AUX/IAA), and glycosyl hydrolase family 9 (Cel3) was associated with lateral root initiation and development 64,65 . The down-regulation of four genes including extensin-2-like and one of chitinases indicated that selenium could affect the growth and development of lateral roots significantly.
Another predominant character is that the degradation of protein occurred in the roots and leaves at the same time but in different tissues and subcellular organelles (Table 4 and Fig. 7G and H). The gene of ubiquitin-conjugating enzyme E2, another key member in the UPP, was also up-regulated in the roots. The gene of V-type H + -transporting ATPase subunit I (VHA-a2), which is located on mature phagosomes, was up-regulated in the leaves. There were more genes responding to the stress from selenium, such as drought tolerance-and pathogen resistance-related genes. These included defensin, pathogenesis-related protein 1, aquaporin PIP, ferritin heavy chain, ARALYDRAFT_483040 (defense response), ATHB-12, transcription factor TGA, ANAC019, RD26, group I late embryogenesis abundant protein (LEA), and protein transport protein SEC. 23 as well as four metal binding proteins, of which only pathogenesis-related protein 1 and aquaporin PIP were down-regulated in the roots and the others were all up-regulated either in the roots or leaves. Surprisingly, the drought tolerance-related genes ATHB-12 48,66 , ANAC019 67,68 , RD26 69 and group I LEA 70 were significantly up-regulated in the leaves, and aquaporin PIP was down-regulated in the roots under the selenium stress. The process of new protein modification was accelerated by up-regulation of Sec. 23, which initiated the COP II coat complex assembly 71 . The pathogen resistance-related genes showed a puzzling change in expression: pathogenesis-related protein 1 was down-regulated in the roots, and transcription factor TGA, whose members interact with the key components (ankyrin repeat protein and non-expresser of pathogen-related (PR) (NPR1) genes) in the SA defense signaling pathway 72 , was up-regulated in the leaves.
Regarding the senescence-associated physiology, the RAV-like factor, which can inhibit the growth of plant leaf, root and stem 73,74 , was up-regulated, and the senescence-associated protein was down-regulated in the leaves.
The SPT2 chromatin protein which was up-regulated in the roots is an important histone chaperone and can facilitate ribosomal DNA transcription through chromatin remodeling 75 . Together with the up-regulation of histone H2A, these results suggest that selenium could function in the process of gene expression in high concentration of Se in the roots. The up-regulation of CCR4-NOT transcription complex subunit 6 also supported this speculation.

Discussion
Selenate is the initial compound of selenium metabolism. Selenium is chemically similar to sulfur and is assimilated by plants via the same metabolic pathways 76,77 . Most plants nonspecifically take up selenate from the environment by means of sulfate transporters and assimilate selenate into organic forms of Se via S metabolic pathways 7 . The conversion of selenate to selenite requires the continuous action of two enzymes (Fig. 1). ATP sulfurylase (APS) mediates the binding of selenate with ATP, forming adenosine phosphoselenate (APSe). This compound is then reduced to selenite through APS reductase (APR) 16 . However, we found that SOX was up-regulated when selenite was added to the culture solution of the C. hupingshanensis seedlings. Therefore, we can deduce that selenite might be converted to selenate first and then was incorporated into ATP by APS, reduced to selenite by APR, and reduced to selenide before finally being incorporated into SeCys in the root tissue of C. hupingshanensis seedlings (Fig. 8).

The storage function of the vacuole plays an important role in selenium tolerance. After
APSe formed, the members of glutathione S-transferase family genes, GST u4 (in leaves) (Figs 8 and 9B,C) were up-regulated to transfer the selenate ion to GSH and to form glutathione-S conjugate (GS-X). ABCC2 (in leaves), belonged to the subfamily C (CFTR/MRP) of ATP-binding cassette superfamily, which proved to be the long-sought and major vacuolar plant PC transporters 78  significantly up-regulated in the leaves (5.26 folds at 100 μg Se/L and 5.19 folds as 80,000 μg Se/L) of C. hupingshanensis seedlings when treated with selenite. Additionally, four genes encoding metal ion binding proteins were up-regulated. Therefore, we deduced that partial selenate was first chelated by glutathione-derived peptides with glutathione sulfur transferase (GST) and was then transported into the vacuole by MRP2 of C. hupingshanensis to detoxify and sequester the heavy metal in the roots and leaves. All these results suggest that the storage of the vacuole is an important way to tolerate selenium in C. hupingshanensis seedlings.
Transamination is an important mechanism of selenium detoxification. The fate of SeCys enormously influences the capacity of selenium tolerance in plants. The misincorporation of SeCys was believed to be the main reason for selenium toxicity 4,11,13,14 . Therefore, the fate of SeCys will determine the toxicity. SeCys can be methylated by SeCys methyltransferase (SMT), or it can be converted into SeMet by CγS, oxidized by SL and specifically broken by CpNifS 16 . The process of SeCys conversion to SeMet still plays an important role, as demonstrated by the up-regulation of CγS, CBL and MET in the leaves. However, the transamination of SeCys or its oxides, such as L-cysteate and L-cysteine-sulfinate, has not been studied in plants. In this study, we found  that the expression of Aat was up-regulated in the roots compared with the control. Furthermore, PLR1 was also up-regulated in the roots. Therefore, we can deduce that SeCys deamination by Aat is an important pathway for detoxification of selenium in the roots of C. hupingshanensis seedlings (Figs 8 and 9A).  Selenation is the more common mechanism for selenium detoxification. Although it is not the only route, selenation, which is similar to sulfation in plants, is another pathway for detoxification of selenium. APR is a key enzyme in both sulfate and selenate reduction 49 which was down-regulated with 3.18-fold in the roots. Therefore, the new metabolic pathway to transfer selenate seemed to be more important. Two ways were found for selenate stress: (1) the selenate is chelated by GSH and then transported into the vacuole; (2) the APSe, 3′-phosphoadenosine 5′-phosphoselenate synthase (PAPSeS) and aryl sulfotransferase (SULT1A), which are present in animals and prokaryotes but were found only in Ricinus communis L. 28,29 , transfer selenate to a phenolic hydroxy group, forming selenocompound substrates 28 . When the concentration of selenium increased to 80,000 μg Se/L, a stress concentration, the gene of SULT1A, whose product is located in the sulfur and glucosinolates metabolism pathway, increased 4.05-fold in the roots (Table 4). This suggests that selenation also occurred, as selenium increased to a high concentration in the roots of C. hupingshanensis seedlings. The response of selenation to selenium is not the same in the leaves. Whether challenged with a low or high concentration of selenium, the expression of SULT1A increased 4.94 or 7.17-fold (Table 2). That is, the response of selenation was more common in leaves when the C. hupingshanensis seedlings were treated with selenium. Therefore, we can conclude that selenation is the more common method of selenium stress tolerance in C. hupingshanensis seedlings (Fig. 8).

Degradation of selenoproteins is important for selenium detoxification. Both of the
Selenomethionine and selenocysteine are seleno-amino acids which can be misincorporated into proteins in plants. Cysteine plays an important role in maintaining the structure and function of proteins, including those involved in catalysis, redox regulation, formation of disulfide bridges and metal binding. The substitution of cysteine with selenocysteine in nonspecific selenoproteins could create either a diselenide bridge or a mixed selenide-sulfide bridge (or selenosulfide bridge) with different properties, and a deformed protein could be formed 13 .
A non-specific selenocysteine incorporated into selenoproteins in other situations and non-specific accumulation of selenomethionine proteins are not considered to be as deleterious as the more reactive selenocysteine, and malformed proteins can be formed 13 . The formation of deformed or malformed selenoproteins induced by chaperone-mediated processes and the proteolysis of irreparable proteins through the lysosome or the ubiquitin-proteasome pathway (UPP) can also occur. The mechanisms of preventing the formation of selenoproteins are related to elevated selenium tolerance in plants 79 . The mechanisms of preventing the formation of selenoproteins are associated with increased selenium tolerance in plants 4 . In this study, the genes of E3 ubiquitin-protein ligase MUL1 and RNF13 (treated with 100 μg Se/L) as well as of ubiquitin-conjugating enzyme E2 7 (treated with 80,000 μg Se/L) were up-regulated with 4.87, 5.55 and 2.82-fold in the roots, respectively. These changes are similar to the observations in Chlamydomonas reinhardtii, Stanleya pinnata and rice 4,79,80 . However, the puzzling change was that the gene of the 26S proteasome regulatory subunit T5 (Rpt5), an important subunit for assembly of the 26S proteasome, was down-regulated with 9.70-fold in all tissues and at both selenium concentrations. On the other hand, regarding phagosome and phagolysosome, the gene of cysteine-type peptidase was up-regulated with 9.98 fold in the roots (treated with 100 μg Se/L), and the gene of VHA-a2 was up-regulated with 10.24 fold in the leaves (treated with 80,000 μg Se/L). All these results indicate that the degradation of selenoproteins plays an important role in selenium detoxification (Figs 8 and 9A,B).

Methods
Plant materials. The seeds of C. hupingshanensis were harvested from the Yutangba Se mining field, which is located on the Enshi area in western Hubei province, China. Plants were grown in a growth chamber with light illumination (ca. 1600 mol −2 ms −1 ) over a 16/8 day and night at 20 ± 2 °C. Hoagland solution was sprayed every two days. Five-month-old C. hupingshanensis seedlings were harvested. The roots were washed first with tap water and then with deionized water (≥18 MΩ, Millipore, Bedford TM , USA) to exclude contamination from the surface. The washed seedlings, which were cultured in 100 mL of Hoagland solution, were divided into three experimental groups based on the results of our previous studies: control, supplementation with 100 μg Se L −1 and supplementation with 80,000 μg Se L −1 , using sodium selenite (analytical reagent, Sinopharm Chemical Reagent Co., Ltd, Shanghai, China). After 24 hours, the C. hupingshanensis seedlings were harvested.

RNA isolation, library preparation and sequencing.
For the transcriptome sequencing, the roots and leaves were separated and analyzed independently. According to the manufacturer's protocol, total RNA was extracted using a Qiagen total RNA isolation system (RNeasy Plant Mini Kit, 74904, Qiagen). The total RNA samples were treated by the following protocol: DNA degraded by DNase I; the oligo (dT) magnetic beads were used for mRNA enrichment; the mRNA was then fragmented into short fragments by mixing with the fragmentation buffer. Then, the cDNA was synthesized using PrimeScript ™ Double Strand cDNA Synthesis Kit (Takara) according to the manufacturer's protocol. The double-strand cDNA was purified with magnetic beads. And then 3′-end single nucleotide A (adenine) addition was performed. Finally, sequencing adaptors were ligated to the fragments. The fragments were enriched by PCR amplification. During the QC step, the Agilent 2100 Bioanalyzer and ABI StepOne Plus Real-Time PCR System were used to qualify and quantify the sample library. The primer information for the real-time quantitative PCR was shown in Table 5. Actin of C. Hupingshanensis (Chp Actin) served as internal controls to normalize the targets for quantification. The libraries were then sequenced on the Illumina HiSeq TM 2000 platform.
Transcriptome assembly and annotation. The raw data were obtained after deep transcriptome sequencing. After quality analysis using the fastqc program (version 0.10.1), the raw data were processed to clip sequencing adapters and filter low-quality reads using Trimmomatic software (version 0.30). The remaining clean reads were used to assemble transcripts using the Trinity program (http://trinityrnaseq.sourceforge.net/) embedded with three individual modules (Inchworm, Chrysalis and Butterfly) which were run consecutively. The parameters for Trinity included the following: -seq-Type fq, -min_contig_length 100, -min_glue 3, -group_ pairs_distance 250, -path_reinforcement_distance 85 and -min_kmer_cov 3. The following parameters were also used in Trinity: min_glue = 1, V = 10, edge-thr = 0.05, min_kmer_cov = 2, path_reinforcement_distance = 150, and group_pairs_distance = 500. The assembled contigs were finally joined together to make scaffolds. The Gene Indices Clustering Tools (TGICL, version 2.1) program was used to form unigenes, and Phrap (http://www.phrap. org/) was used to assemble the scaffolds and cluster them.
To annotate these unigenes, all sequences were subjected to blastx alignment (e-value < 1e-5) with the NR (non-redundant protein sequence database, release 20130408), Swiss-Prot (release 2013_03), KEGG (Kyoto Encyclopedia of Genes and Genomes) and COG (Clusters of Orthologous Groups) databases. The results from these alignments were used to determine the direction of these sequences. Lastly, the sequences that were not aligned to any database were subjected to a ESTScan analysis. For function annotation, a local Blast was used to search against the NT (NCBI nucleotide database), NR, Swiss-Prot, KEGG and COG databases. Blast hits from the NR database was then used to determine gene ontology (GO) terms of unigenes. The Blast2GO (http://www. blast2go.com/b2ghome) program was run to obtain the GO terms.
Differential gene expression analysis. Expression profiling of unigenes was performed. The gene expression levels were quantified using the software package RNASeq by Expectation Maximization (RSEM). The fragments per kilobase per million reads (FPKM) method, which is able to eliminate the influence of different gene lengths and sequencing discrepancy on the calculation of gene expression, was used to estimate the expression level of each gene. Therefore, the gene expression levels normalized as FPKM values can be directly used to assess the differences in gene expression among samples. The FPKM value was calculated according to the length of each unigene and the number of reads mapped to the gene. Statistical analysis was performed to identify differentially expressed genes (DEGs). The false discovery rate (FDR) was calculated to adjust the p-value threshold in the expression analysis. If the FDR was small and the fold change was large, the difference in expression between the two samples was large. The criteria used to detect DEGs were FDR ≤ 0.001 and fold change (Se treated/not treated) ≥1 or ≤−1. In addition, the GO and KEGG pathway analyses were performed for the DEGs similar to the method described above.