Seeds used in the experiments were kindly provided by Genebank Gatersleben of the Leibniz Institute of Plant Genetics and Crop Plant Research. Barley seeds (Hordeum vulgare cv. Morex) were germinated at 20°C in the dark on filter paper moistened with deionized water. After 48 h the primary root length was measured for seeds that were then transferred to polypropylene pots containing 1L of test medium (changed every 24 h to maintain composition). Test media were aerated throughout the exposure time. Relative humidity was 50%, and day/night temperatures were 20/16°C during 16/8 h photoperiod.

The basal nutrient solution composition was KH₂PO₄, 0.4 mM; K₂SO₄, 0.4 mM; MgSO₄ × 7H₂O, 0.6 mM; NH₄NO₃, 1 mM; Ca(NO₃)₂ × 4H₂O, 2 mM; FeNaEDTA, 75 μM; MnCl₂ × 4H₂O, 7 μM; ZnCl₂, 3 μM; CuSO₄ × 5H₂O, 0.8 μM; H₃BO₃, 1.6 μM, and Na₂MoO₄ × 2H₂O, 0.83 μM. The pH was adjusted with 1 M KOH to 5.5–6.0 (Podar, 2013). Test media, except control, were prepared by adding different volumes of stock solutions of CdCl₂ to adjust the concentration.

To determine Cd accumulation in barley cv. Morex, plants were cultivated in control medium and medium supplemented with 80 μM cadmium, a previously determined effective concentration (EC₅₀) causing a 50% inhibitory effect for root growth (Kintlová et al., 2017). Samples were collected and root length measured before and after 5 days of treatment. Both root and shoot tissue samples were washed first in 0.5 M EDTA to remove unabsorbed cadmium, then washed in distilled water. They were further dried in a hot air oven (50°C, 72 h). The cadmium concentration was then measured by atomic absorption spectrometry (AAS).

For transcriptome analysis we used the dataset available in NCBI SRA (sequence read archive) under accession numbers SRR5452097 to SRR5452108. These data were obtained from the barley plants grown in media supplemented with cadmium concentration that corresponds to established EC₅₀ values as well as control plants for each treatment (Kintlová et al., 2017). Sequenced Read quantification of each dataset was performed using Salmon (Patro et al., 2017) and 63,658 H. vulgare CDS (Mascher, 2019) and 3,349,186 repeats (IBSC, 2016). Differential expression analysis was then conducted using R and the package DESeq2 (Love et al., 2014). Gene ontologies were retrieved from the available barley gene annotation (Mascher, 2019) and were further analyzed with the DE genes. Bioinformatic analysis of RNA-seq data was carried out in the computing and storage facilities MetaCentrum infrastructure.

The mRNA levels for selected differentially expressed transcripts was validated by qRT-PCR. As representatives, most up- and down-regulated transcripts from both root and shoot tissues were selected for analysis (Supplementary Table 3). Plant material was grown and treated as described above. The RNA from control and Cd-treated plants (80 μM) was isolated in triplicate by NucleoSpin RNA Plant kit (Macherey-Nagel, Germany) and transcribed to cDNA by Transcriptor high fidelity cDNA synthesis kit (Roche). The qRT-PCR reactions were performed on a Light Cycler 96 real-time PCR system (Roche, USA) using the 2x SYBR Master Mix (Top-Bio, Czech Republic). Relative expression of shoot and root transcripts was assessed using CFX Manager™ Software (Bio-Rad, Hercules, California, USA).

Amino acid identity for each HvPCR2 gene and previously described Arabidopsis PCR genes (Song et al., 2004, 2010) was determined using alignment algorithm MAFFT v7.450 (Katoh and Standley, 2013). The genomic sequence of Hordeum vulgare (IBSC_v2), Oryza sativa Indica group (ASM465v1) Zea mays (B73_RefGen_v4) and Musa acuminata (ASM31385v1) were used to analyze genomic loci containing PCR2 homologs, data were retrieved and visualized using EnsemblPlants website (Howe et al., 2020).

To test the function of HvPCR2 genes, all five identified copies were expressed in Saccharomyces cerevisiae cadmium hypersensitive mutant strain DTY168 (MATα, his6, leu2-3−112, ura3-52, ycf1::hisG) lacking the yeast cadmium factor 1 gene (Szczypka et al., 1994). HvPCR2.1-5 coding sequences were amplified from Morex cDNA (Supplementary Table 4) and cloned into p426 under the control of yeast GDP (glyceraldehyde-3-phosphate dehydrogenase) promoter (Mumberg et al., 1995). Vectors containing HvPCR2 fragments and empty vector p426 were introduced into the DTY168 strain by lithium acetate. Yeast transformants were grown overnight, OD600 was adjusted to 0.5, four 10-fold dilutions were made, and each culture was spotted on SD-Ura plates and SD-Ura plates supplemented with 50 μM CdCl₂. Plates were incubated for 5–6 days at 28°C.

The effect of 80 μM Cd on barley after a 5-day treatment comprised of root shortening and overall growth retardation (Figure 1A). Root growth of Cd-treated plants was inhibited by ~50% in comparison with non-treated plants (Figure 1B; Supplementary Table 1). Further, we measured the Cd content in control and treated plants in order to obtain information about the distribution of Cd in barley plants in both root and shoot tissues. The Cd content in roots was in control conditions approximately double of the Cd concentration in shoots. Cd-treated plants accumulated roughly 6,600 times more Cd in roots than non-treated plants. This result shows a strong accumulation of Cd in the root and restriction of Cd translocation to the aerial tissues in high Cd stress (Figures 1C,D; Supplementary Table 2).

The labeled thymidine analog (EdU) was used to monitor cell cycle abnormalities of Cd treated plants and to determine the Cd effect on root growth at the cellular level. In control plants, we observed an approximately equal number of cells before, during and after the S phase (Figure 2A). Cd-treated roots contained a significantly lower number of S phase cells (Figure 2B). Such a decrease may suggest that the cells entering the S phase are the most susceptible to Cd treatment. After removing the Cd by washing the roots, the number of cells in each state went back to the level similar to the non-treated control (Figure 2C). These findings indicate that the major reduction in root growth is caused by the suppression of the cell division process namely in G1/S cells.

The entire dataset subjected to differential expression analysis contains 20.7 million of non-treated and 24.3 million of cadmium treated RNA reads (Kintlová et al., 2017). The analysis was performed using the DESeq2 package from Bioconductor (Love et al., 2014). We controlled the clustering of both control and cadmium replicates using multidimensional scaling for each tissue (Supplementary Figure 1). We identified 14,459 differentially expressed (DE) transcribed sequences in both root and shoot tissues using a p-value threshold of 0.05 (Figure 3). The overall count of DE after cadmium treatment is 14,459 (Figure 4; Supplementary Tables 6, 7) with 2,927 transcripts that have a significant change in their expression in both tissues while 7,355 were differentially expressed only in the root and 4,177 in the shoot. Separating the set of DE transcripts between up- and down-regulated reveals that most of the transcripts in the root and shoot are down-regulated (3,509 up vs. 3,846 down in root and 1,825 up vs. 2,352 down in shoot). Differentially expressed transcripts present in both tissues follow the same pattern of expression change with more down-regulated transcripts than up-regulated ones. From the 2,927 DE transcripts in both tissues, 903 are up-regulated, and 1,438 are down-regulated in roots and shoots. Nevertheless, 288 transcripts are down-regulated in roots while up-regulated in shoots, and 298 transcripts are up-regulated in roots while down-regulated in shoots (Figure 4).

Analysis of the transcripts related to repeated elements revealed that 2,285 elements (from total 14,459 DE transcribed sequences) are differentially expressed (Supplementary Figure 2). 59% are specifically differentially expressed in the root while 32% in the shoot. About 75% of the total of DE transposable elements are long terminal repeats (LTR). Transposable element related transcripts are more expressed in roots than shoots, similar to the genes-related ones, but unlike gene transcripts, they are more often up-regulated. Most of the expressed LTR belong to an unknown subfamily; otherwise, DE LTR Gypsy are more abundant than LTR Copia, followed by CACTA, nonLTR/LINE, and MITE elements (Supplementary Figure 3; Supplementary Table 6).

The 10 most abundant gene ontology (GO) annotations for both roots and shoots were assessed for DE transcripts. Comparison of the results reveals that regardless of the tissue, most of the genes influenced by cadmium are: involved in oxidation-reduction processes and phosphorylation, have mainly ATP binding or transferase activity functions and are part of the membrane. Moreover, the genes involved in protein phosphorylation and protein kinase activity are also influenced, but the majority of them are up-regulated while the genes involved in catalytic activity are mainly down-regulated (Figure 5).

In order to validate RNA-seq results, we selected a total of six transcripts for which we designed qPCR primers (Supplementary Table 3). HORVU.MOREX.r2.7HG0558200 is up-regulated in shoots, HORVU.MOREX.r2.6HG0456470 is down-regulated in shoots and HORVU.MOREX.r2.1HG0006620 is not differentially expressed in shoots. HORVU.MOREX.r2.2HG0162330 is up-regulated in roots, HORVU.MOREX.r2.7HG0538720 is down-regulated in roots and HORVU.MOREX.r2.7HG0561140 is not differentially expressed in roots (Supplementary Figure 4). The qPCR results correspond to the mRNA level identified by differential expression analysis for each of the six candidates (Supplementary Figure 5).

Among the most Cd upregulated transcripts, we identified homologs of PLANT CADMIUM RESISTANCE 2 (PCR2) gene. Predicted protein sequences show the highest similarity with AtPCR2 protein (Supplementary Table 5), therefore we have designated these genes HvPCR2.1 (HORVU.MOREX.r2.1HG0071220), HvPCR2.2 (HORVU.MOREX.r2.5HG0439660), HvPCR2.3 (HORVU.MOREX.r2.5HG0439680), and HvPCR2.4 (HORVU.MOREX.r2.5HG0439700). While the mRNA level of HvPCR2.1 is significantly upregulated by Cd in both root and shoot tissue, HvPCR2.2-4 are upregulated only in the shoot (Figure 6A). Surprisingly, HvPCR2.2-4 were found to be located in the same locus on barley chromosome 5H (Figure 6B). The fifth copy of HvPCR2 gene (HORVU.MOREX.r2.5HG0439690; denominated HvPCR2.5) is also physically present in this ~35 kb spanning region as well (Figure 6B). However, HvPCR2.5 does not exhibit significant upregulation upon Cd treatment. Similar tandem arrangement of PCR2 homologs have been found in rice (~30 kb locus on chromosome 3 containing six PCR2 homologs), corn (~90 kb locus on chromosome 5 containing three PCR2 homologs) and banana (~30 kb locus on chromosome 3 containing three PCR2 homologs) (Supplementary Figure 23).

To elucidate the function of HvPCR2 genes in metal tolerance, we expressed them in the cadmium sensitive yeast strain DTY168 (Δycf1), which showed restricted growth at increased Cd concentrations. Growth of the DTY168 harboring empty vector p426 GDP was strongly inhibited by 50 μM CdCl₂, whereas the same yeast mutant strain expressing either of the HvPCR2 genes was able to restore Cd tolerance (Figure 6C).

To provide quantitative information about the presence of cadmium in barley plants, we measured the Cd concentration in roots and shoots separately (Figures 1C,D). We found that under cadmium treatment roots and shoots were able to accumulate a high concentration of this element. While Cd-treated shoots accumulate up to 300 times more Cd than in the non-treated condition, roots accumulated more than 6,600 times of this element. This difference indicates the limitations in the transport of cadmium in the aerial part of the plant.

Monitoring the cell cycle by EdU labeling revealed a lower number of S phase cells in Cd-treated roots. Such a decrease may suggest that the cells before G1/S transition are the most susceptible to Cd treatment. The recovery of the number of cells after removing the Cd by washing the root indicates that the major reduction in root growth is caused by the suppression of the cell division process namely in G1 phase cells (Figure 2). While the negative influence of Cd toxicity on cell proliferation is well-known in barley and wheat (Zhang and Yang, 1994; Paradiso et al., 2008), this phenomenon is not yet fully understood. Presumably, a Cd-induced ROS accumulation may cause an oxidative posttranslational modification of cyclin D and A-type cyclin-dependent kinase, thereby reducing the functionality of essential G1/S transition regulators (Pena et al., 2012). Similar reductions in the number of S phase cells have been reported in Cd-treated Zea mays (Bertels et al., 2020), Sorghum bicolor (Zhan et al., 2017), and Vicia faba (Zabka et al., 2021) as well as in Al-treated barley (Jaskowiak et al., 2018). Although Cd is also known to cause abnormal mitosis and aberrant chromosomes in barley (Shi et al., 2016), our data suggest a substantial effect of Cd toxicity at the level of G1-to-S-phase transition.

In our study, we identified over 14,400 transcripts differentially expressed in Cd treated plants, with about half of them being down-regulated. The total number of DE transcripts were lower in shoots (7,104) than in roots (10,282), consistent with the Cd concentration level detected in both tissues (Figures 1C,D). Previous transcriptomic studies reporting the effect of cadmium on barley focused on Cd tolerance, rather than response to Cd excess (Cao et al., 2014; Derakhshani et al., 2020). As the general objective of our study was to describe the overall response to cadmium in barley, we identified higher number of differentially expressed transcripts.

First, we analyzed the repetitive fraction of the barley genome at the transcriptional level to quantify the impact of cadmium on the activity of individual elements. Abiotic stress is frequently accompanied by the transcriptional and transpositional activation of transposable elements (TEs) (Grandbastien, 1998). Moreover, cadmium related stress leads to the establishment of a new balance of expressed/repressed chromatin (Greco et al., 2012). We found that the vast majority of transcriptionally affected repeats represent LTR retrotransposons (Copia, Gypsy and unknown families) and CACTA, nonLTR/LINE and MITE elements (Supplementary Figure 3). Transcriptional response to cadmium treatment is more evident in root tissue (Supplementary Figures 2, 3) including both upregulation and downregulation of individual TEs (Supplementary Table 6).

DE transcripts are mostly connected with oxidation-reduction processes, protein phosphorylation, metabolic processes, regulation of transcription and transmembrane transport (Figure 5). Among them, the most upregulated transcripts match genes encoding proteins involved in ROS metabolism and oxidative stress such as tau class glutathione-S-transferase (GST, Supplementary Figures 6, 7), purple acid phosphatase (PAP, Supplementary Figure 8), multidrug resistance-associated protein 3 (MRP3, Supplementary Figure 9), and cytochrome P450 enzymes (P450, Supplementary Figures 10, 11). All of these are also known to play a role in the detoxification of xenobiotics or abiotic stress signaling (Bovet et al., 2003; Schenk et al., 2013; Kumar and Trivedi, 2018; Pandian et al., 2020). Cadmium is known to induce oxidative stress in plants (Grobelak et al., 2019); therefore, the enhancement of ROS scavenging mechanisms represents a straightforward reaction to elevated Cd concentration. In particular, increased expression of GST and P450 enzymes in response to Cd stress has been previously described in rice (Ogawa et al., 2009) and the simultaneous overexpression of both enzymes in Medicago truncatula led to improved tolerance to mercury contamination (Zhang et al., 2013). In addition, four genes encoding homologs to putative chromatin remodeler OXIDATIVE STRESS 3 (OX3) were upregulated (Supplementary Figure 12). Although the molecular mechanism of OX3 action is not completely understood, this enzyme provides tolerance to Cd and oxidative stress in Arabidopsis and fission yeast (Blanvillain et al., 2009).

Other coding genes showing significant upregulation after Cd treatment, include those involved in cell wall formation and modification. Transcripts coding for lignin biosynthesis components, namely phenylalanine ammonia lyase (PAL; Supplementary Figure 13), shikimate O-hydroxycinnamoyltransferase (HCT; Supplementary Figure 14), cinnamoyl-CoA reductase (CCR; Supplementary Figure 15), cinnamyl alcohol dehydrogenase (CAD; Supplementary Figure 16), laccase 7 (LAC7; Supplementary Figure 17), and peroxidase (POD; Supplementary Figure 18) have been highly upregulated in Cd response. PAL is the key enzyme of the lignin biosynthesis acting at the beginning of this metabolic pathway. Accordingly, the Arabidopsis pal1 pal2 pal3 pal4 quadruple mutant showed approximately four times decreased lignin content in its cell wall compared with the wild-type (Huang et al., 2010). PAL activity or expression is known to be induced by Cd stress in various plant species including pea (Głowacka et al., 2019), Medicago sativa (Gutsch et al., 2019), Salix matsudana (Yang et al., 2015), rice (Hsu and Kao, 2005), and wheat (Shakirova et al., 2016). HCT, CCR and CAD catalyze the synthesis of monolignols, which are then polymerized into lignin by LAC and POD enzymes (Liu et al., 2018). As the lignification is improving the cell wall impermeability and lignin itself binds cadmium (Parrotta et al., 2015), the enhanced expression of these enzymes is likely to result in a lower Cd uptake and a higher Cd sequestration capacity. Notably, altered lignin biosynthesis was proposed to contribute to cadmium tolerance in Zn/Cd hyperaccumulator Noccaea caerulescens (van de Mortel et al., 2006, 2008). Other highly Cd responsive transcripts encode homologs of wound induced protein 1 (HORVU.MOREX.r2.1HG0053620; Supplementary Figure 19) and Leucine-rich repeat receptor-like protein kinase (HORVU.MOREX.r2.2HG0164280, Supplementary Figure 19), both of which have been reported to function in cell wall integrity sensing and responding to biotic and abiotic stress (Logemann and Schell, 1989; Van der Does et al., 2017).

Metal ion transporters play significant role in heavy metal homeostasis and detoxification (Jain et al., 2018). In our dataset, five deregulated transcripts have been identified to encode the members of the ZRT/IRT-like protein (ZIP) family (Supplementary Figure 20). ZIP transmembrane proteins are known to mediate cellular uptake, intracellular trafficking and detoxification of heavy metal cations such as Zn, Fe, Cd, Cu, Mn, Co, and Ni (Ajeesh Krishna et al., 2020). The role of selected ZIP transporters has been studied extensively in monocots (Li et al., 2013; Evens et al., 2017; Zheng et al., 2018) including barley (Pedas et al., 2009; Tiong et al., 2015). Interestingly, QTL mapping in barley identified two ZIP transcripts among the candidate genes for Cd tolerance (Derakhshani et al., 2020). In our dataset, we identified three previously studied ZIP transporters–HvZIP1 (HORVU.MOREX.r2.3HG0273580), HvZIP3 (HORVU.MOREX.r2.2HG0158440), and HvZIP6 (HORVU.MOREX.r2.1HG0020460). HvZIP1 was found to be upregulated by Cd in both tissues, but the effect is much stronger in roots (Supplementary Figure 20). An identical expression pattern has been found for OsZIP1, a rice ortholog for HvZIP1, that functions as a heavy metal efflux transporter maintaining detoxification of Zn, Cu and Cd in rice (Liu et al., 2019). Upon cadmium exposure, the HvZIP3 shows upregulation in roots but downregulation in shoots (Supplementary Figure 20). The transgenic barley plants with reduced HvZIP3 expression showed a significant increase in grain Cd accumulation (Sun et al., 2015). Further, a co-expression of rice HvZIP3 homolog (OsZIP3) with rice heavy metal transporter genes HEAVY METAL ATPASE 2 (OsHMA2) and LOW-AFFINITY CATION TRANSPORTER1 (OsLCT1) resulted in a decrease in Cd accumulation and enhanced Cd tolerance (Tian et al., 2019). HvZIP6 shows strong upregulation in shoots only (Supplementary Figure 20). So far, HvZIP6 was studied only in relation to Zn homeostasis (Tiong et al., 2015) but OsZIP6 is known to transport Cd²⁺ cations (Kavitha et al., 2015), therefore HvZIP6 may take part in shoot Cd detoxification. Two previously unknown barley ZIP transporters encoded by HORVU.MOREX.r2.2HG0098160 and HORVU.MOREX.r2.7HG0603680, show downregulation after Cd treatment (Supplementary Figure 20), which may suggest these could mediate Cd uptake.

Another group of deregulated genes encoding cation transmembrane transporters are members of the NATURAL RESISTANCE ASSOCIATED MACROPHAGE PROTEIN (NRAMP) family. In Cd-treated roots, HvNRAMP5 (HORVU.MOREX.r2.4HG0337880) exhibited several-fold downregulation (Supplementary Figure 21), which is consistent with the role of HvNRAMP5 as a Cd uptake transporter (Astolfi et al., 2014; Wu et al., 2016). Conversely, the HvNRAMP5 mRNA level slightly increased after Cd excess in shoot, which may imply involvement of this gene in shoot Cd distribution in stress condition. The HvNRAMP2 (HORVU.MOREX.r2.4HG0331160) were downregulated in shoots (Supplementary Figure 21). As the genome-wide association mapping study has proposed OsNRAMP2 to be responsible for Cd accumulation (Zhao et al., 2018), the downregulation of this transcript in barley might serve as a defense against accumulation of high Cd concentrations in shoot. The two remainder transcriptionally deregulated transporters from the NRAMP family were HvNRAMP1 (HORVU.MOREX.r2.7HG0610240; upregulated in roots) and HvNRAMP6 (HORVU.MOREX.r2.3HG0215920; upregulated in shoots) (Supplementary Figure 21).

Two transcripts encoding HEAVY METAL ATPASE (HMA), a P-type ATPase cation transporter family, were identified as being responsive to Cd. HvHMA1 (HORVU.MOREX.r2.5HG0401460) is expressed only in shoot tissue and no transcripts are detected after Cd treatment (Supplementary Figure 22). This pattern is similar to results published previously (Mikkelsen et al., 2012) and suggests that Cd induced downregulation of protoplast- and aleurone layer-localized HvHMA1 helps to prevent the cadmium transport into grains and protoplasts. The mRNA level of HvHMA3 gene (HORVU.MOREX.r2.7HG0603650) was decreased by Cd treatment (Supplementary Figure 22). Unfortunately, the function of this HMA3 copy has not been studied yet, so it is hard to clarify the effect of its downregulation. On the other hand, different HvHMA3 gene (HORVU5Hr1G094430) has been proposed to play a crucial role in grain Cd accumulation (Wu et al., 2015; Lei et al., 2020).

Remarkably, four homologs of PLANT CADMIUM RESISTANCE 2 (PCR2), encoding a Cys-rich domain containing membrane Cd/Zn efflux transporter (Song et al., 2010), have been identified among the Cd upregulated transcripts in shoot (Figure 6A). PCR1 and PCR2 genes were previously studied in Arabidopsis and it was proposed that they play a role in Zn/Cd detoxification and distribution (Song et al., 2004, 2010). The Arabidopsis pcr2 loss-of-function mutant is sensitive to Cd (Song et al., 2010). In rice, overexpression of OsPCR1 and OsPCR3 mediated Cd tolerance and lowered Cd accumulation (Wang et al., 2019). As there are no data regarding the function of PCR2 genes in monocots, we expressed barley HvPCR2 genes in a Cd sensitive yeast mutant. All five HvPCR2 genes restored the Cd tolerance in yeast (Figure 6C), hence it is probable that they retain a similar function as Arabidopsis homologs. Moreover, we have found out that four HvPCR2 genes are located on barley chromosome 5H in tandem arrangement (Figure 6B), probably due to local gene duplication events. A similar arrangement of multiple PCR2 copies was found in genomic regions of rice, corn and banana (Supplementary Figure 23), suggesting that PCR2 multiplication might be common in monocots. Duplication of heavy metal transporter genes is known to contribute to Cd tolerance in metal hyperaccumulators Arabidopsis halleri and Noccaea caerulescens (Hanikenne et al., 2008; Ó Lochlainn et al., 2011). Beside increasing gene dosage, the duplication may also allow functional divergence. It is therefore possible that the increased copy number of PCR2 genes originated as an adaptive response to Cd toxicity and/or Zn homeostasis requirements.