Global leaf and root transcriptome in response to cadmium reveals tolerance mechanisms in Arundo donax L

The expected increase of sustainable energy demand has shifted the attention towards bioenergy crops. Due to their know tolerance against abiotic stress and relatively low nutritional requirements, they have been proposed as election crops to be cultivated in marginal lands without disturbing the part of lands employed for agricultural purposes. Arundo donax L. is a promising bioenergy crop whose behaviour under water and salt stress has been recently studied at transcriptomic levels. As the anthropogenic activities produced in the last years a worrying increase of cadmium contamination worldwide, the aim of our work was to decipher the global transcriptomic response of A. donax leaf and root in the perspective of its cultivation in contaminated soil. In our study, RNA-seq libraries yielded a total of 416 million clean reads and 10.4 Gb per sample. De novo assembly of clean reads resulted in 378,521 transcripts and 126,668 unigenes with N50 length of 1812 bp and 1555 bp, respectively. Differential gene expression analysis revealed 5,303 deregulated transcripts (3,206 up- and 2,097 down regulated) specifically observed in the Cd-treated roots compared to Cd-treated leaves. Among them, we identified genes related to “Protein biosynthesis”, “Phytohormone action”, “Nutrient uptake”, “Cell wall organisation”, “Polyamine metabolism”, “Reactive oxygen species metabolism” and “Ion membrane transport”. Globally, our results indicate that ethylene biosynthesis and the downstream signal cascade are strongly induced by cadmium stress. In accordance to ethylene role in the interaction with the ROS generation and scavenging machinery, the transcription of several genes (NADPH oxidase 1, superoxide dismutase, ascorbate peroxidase, different glutathione S-transferases and catalase) devoted to cope the oxidative stress is strongly activated. Several small signal peptides belonging to ROTUNDIFOLIA, CLAVATA3, and C-TERMINALLY ENCODED PEPTIDE 1 (CEP) are also among the up-regulated genes in Cd-treated roots functioning as messenger molecules from root to shoot in order to communicate the stressful status to the upper part of the plants. Finally, the main finding of our work is that genes involved in cell wall remodelling and lignification are decisively up-regulated in giant reed roots. This probably represents a mechanism to avoid cadmium uptake which strongly supports the possibility to cultivate giant cane in contaminated soils in the perspective to reserve agricultural soil for food and feed crops.


Introduction
Nowadays, globe climate change has become one of the most urgent problems humans have to deal with due to the ongoing Green House Gas (GHG) emissions into the atmosphere. Fossil fuels based on organic origin and non-renewable energy sources are by now unsustainable, thereby new energy policies should be considered in order to meet the global energy demand [1]. Among the different renewable energy sources, second generation of biofuels (biomass) is becoming an interesting candidate because of its peculiarity to be converted in a variety of products, such as solid fuels, liquid fuels, heat, electricity, and hydrogen. In particular, biomass is mainly produced from organic origin by using perennial herbaceous plants and fast-growing trees, which are able to limit GHG levels [2,3]. Bioenergy crops have the capability to be cultivated on marginal land, this fact being very attractive in order not to compete with feed and food crops for land use [4]. This fact assumes a crucial importance in the perspective that the global world's population is expected to reach 9.7 billion by 2050 [5], thus determining a considerable increase for agricultural land demand [6]. Such marginal sites have little or no agricultural or industrial value, are characterized by little potential for profit, and often have poor soil or other undesirable characteristics such as high salt content, inadequate water supply or heavy metal (HM) contaminations. Although some HMs are categorized as essential elements because of their positive impact on the plant growth and crop yield (e.g., B, Cu, Co, Fe, Mn, Ni, Zn), some of them, such as Cd, Hg, Pb and As, do not play any role in plant metabolism and can reduce the crop productivity in case their levels reach toxic concentrations [7]. Given both its high toxicity and the high capacity of diffusion in lands, cadmium has led to contaminated soils worldwide [8,9] and in particular in Europe where its concentration in the topsoil raised from < 0.01 to 14.1 ppm in the last years [10]. Moreover, due to its high solubility in soil, cadmium can be easily extracted by the plant root system [11], from where it is accumulated in the above-ground regions inhibiting plant growth and causing a threat to animal and human health through the food chain [12]. Cadmium accumulation-related molecules include many transporters, chelators, some amino acids and organic acids, and the genes encoding the corresponding proteins/enzymes were functionally characterized in both A. thaliana and O. sativa [13]. A group of transmembrane proteins involved in metal transport and homeostasis is represented by the family of NRAMP metal ion transporters which is supposed to be the major transporter family of Cd 2+ from the soil into root cells [14]. Once inside the cells, the mechanisms by which Cd exerts toxicity include: a) substitution of some essential metal ion (e.g. Zn 2+ and Fe 2+ ) or blocking functional groups which leads to inactivation of biomolecules [15], b) a tight binding with thiol groups of proteins which destroys their structure and function [16], and c) generation of reactive oxygen species (ROS), which brings to oxidative stress [17,18].
Arundo donax L. also known as giant reed, is a perennial rhizomatous grass species, genus Arundo, belonging to Poaceae family. Among rhizomatous grasses dedicated to energy production, A. donax represents one of the most promising bioenergy crops [19,20] because of its high biomass production, both low irrigation and nitrogen input requirements, and its high tolerance to abiotic stress conditions, including herbicide, salinity and heavy metals [21,22]. Notably, it has been proposed as species to be employed for phytoremediation [23] due to its ability to accumulate and tolerate high doses of heavy metals, such as Ni, Cd and As [24,25].
Transcriptomic analysis represents a powerful tool to elucidate the molecular mechanisms by which plants accumulate, translocate and detoxify HMs [26]. Moreover, RNA Sequencing (RNA-Seq) and de novo assembly of the transcriptome allow the discover and the quantitative determination of all the expressed genes for those species whose genome sequence is not available yet, such as A. donax. Existing genomic resources of A. donax were provided by RNA sequencing extracted from different organs (leaf, culm, bud and root) of giant reed ecotypes subjected to either normal [27][28][29] or under water stress condition [30]. More recently, the analysis of the transcriptional response of giant reed was performed after a long-term period of salt stress in two different A. donax ecotypes [31,32]. Overall, these studies contributed both to the drafting of candidate gene list that can be used in molecular breeding projects and in revealing the role of plant genes in the abiotic stress response. However, the transcriptional response of A. donax subjected to cadmium treatments is still incomplete and the molecular mechanism of cadmium stress on giant reed metabolic processes remains poorly understood. Recently, Shaheen et al. [33] evaluated the effects of increasing concentrations of cadmium on the expression of selected genes (carotenoid hydroxylase, amidase, glutathione reductase, bHLH, NRAMP and YSL) in Arundo donax L. cultivated in hydroponic solution. The highest expression for these genes was observed in plants exposed to the highest Cd concentration (100 mg/L). Moreover, the activity of several enzymes involved in the ROS scavenging (SOD, CAT, POD) was also measured revealing their activation at the highest cadmium concentration and confirming the onset of a secondary oxidative stress [33]. A study conducted in cohorts of A. donax strongly suggests that phytochelatin encoding genes (AdPCS1, AdPCS2, and AdPCS3) most likely contribute to Cd detoxification in A. donax and the presence of multiple PCS isoforms seems to be advantageous both to provide higher general levels of phytochelatin biosynthesis and to increase flexibility in HM resistance [34].
In this study, considering the frequency of cadmium in soils of Mediterranean basin and the perspective of potential use of marginal land to cultivate bioenergy crops, we sequenced and de novo assembled the giant reed leaf and root transcriptomes after a prolonged period of cadmium treatment by using RNA-Seq technique. The aim of our work was to gain novel insight into the cadmium stress tolerance and to shed new light on the distinct role of leaf and root in the dynamics of heavy metal stress response in plants.

Plant material and application of cadmium nitrate
The experiment was conducted at the Department of Agricultural, Food and Environment (Di3A) of the University of Catania, using G10 ecotype coming from Fondachello (Italy) (N 37 o 45, E 15 o 11'), collected for the Giant reed Network project [35]. The trial started on May 4 th 2020, by filling each pot with 8 kg of clay soil and kept in open air. The contamination of each pots was achieved by using 4 ppm (4 mg/kg) of cadmium nitrate Cd(NO 3 ) 2 for treated plants and 5.774 g of NH 4 NO 3 for control plants, so as to equilibrate the N concentration between the two conditions. The aforementioned concentration is higher than the reference threshold established by the European Union directive 86/278/EEC on Environment protection which establishes the maximum concentration of cadmium on soil should range between 1 and 3 ppm [36]. Cadmium contents above 3 mg/kg are generally thought to indicate contaminated soil [37]. Afterwards, one litre of tap water was added to each pot to allow the element adsorption by soil colloids. Successively, in each pot a single A. donax stem node was transplanted and the irrigation was performed three times a week by adding one litre of tap water to avoid Cd(NO 3 ) 2 leaching. The pots were arranged according to a randomized block factor scheme, considering three biological replicates for each treatment, with a total of six experimental units. Before sampling (July 28 th , 2020, after three months from the start of the trial), the following morphobiometric and physiological parameters were evaluated: number of culms per pots, height of the main culm, net photosynthesis efficiency measured by LCi-T portable photosynthesis system (ADC BioScientific Ltd). Moreover, the measurement of biomass was performed [35]. R software was used for standard deviation analysis. The data were submitted to one-way ANOVA test using R Studio. The averages with pvalue ≤ 0.05 were separated by the TukeyHSD test.

Sample collection and RNA extraction
On July 28 th 2020, fully expanded, no senescing G10 leaves and roots were harvested and immediately frozen with liquid nitrogen [32]. RNA isolation was carried out by using the Spectrum Plant Total RNA Extraction kit (Sigma-Aldrich, St. Louis, MO, USA) according to the manufacturer's instructions. RNA degradation and contamination were monitored by electrophoresis with 1% agarose gel. RNA purity and concentration were assayed using the NanoDrop spectrophotometer (Ther-moFisher Scientific, Waltham, MA, USA) [32]. Before to be sequenced, the RNA samples were subjected to quality parameter evaluation. RNA integrity was assessed using the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA) [31].

Library preparation for transcriptome sequencing
One µg of RNA was used as input material for library preparations (twelve libraries, one for each sample). Sequencing libraries were generated using NEBNext ® Ultra ™ RNA Library Prep Kit for Illumina ® (New England Biolabs, Ipswich, MA, USA) following manufacturer's recommendations [31]. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer (5X). First strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase (RNase H) as synthesizing enzyme [31]. Second strand cDNA synthesis was subsequently performed using RNase H to insert breaks into the RNA molecule and DNA Polymerase I as synthesizing enzyme. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of 3′ ends of DNA fragments, NEBNext Adaptor with hairpin loop structure were ligated to prepare for hybridization. In order to select cDNA fragments of preferentially 150 ~ 200 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, MA, USA). Then 3 μl USER Enzyme by NEB were used with size-selected, adaptor-ligated cDNA at 37 °C for 15 min followed by 5 min at 95 °C before PCR. Then PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer. Finally, PCR products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system [31].

Clustering and next generation RNA sequencing
Cluster generation and sequencing were performed by Novogene (UK) company Limited (25 Cambridge Park, Milton Road, Cambridge, CB4 OFW, United Kingdom). The clustering of the index-coded samples was performed on a cBot Cluster Generation System using a PE Cluster kit cBot-HS (Illumina) [31]. After cluster generation, the library preparations (twelve libraries, one for each sample) were sequenced on Illumina HiSeq2000 platform to generate paired-end reads whose size was paired-end 2 × 150 bp reads. Raw reads in fastq format were firstly processed through in-house perl scripts. In this step, clean data were obtained by removing reads containing adapters, reads containing poly-N and lowquality reads [31]. At the same time, Q20, Q30, GC-content and sequence duplication level of the clean data were calculated. All the downstream analyses were based on clean data with high quality (Table 1) [31].

Quantification of gene expression and differential expression analysis
Gene expression level was estimated by RSEM software (1.2.28 version) by mapping back each clean read onto assembled transcriptome and read counts for each gene were then obtained from the mapping results [47]. Furthermore, the read counts of each gene have been used as input data for DESeq2 (1.26 version, padj ≤ 0.05), to obtain differentially expressed genes (DEGs) [48]. Comparisons were made to identify the set of differentially expressed genes between control (CK) and cadmium (Cd) treatments for each tissue (Cd_L_vs CK_L; Cd_R vs CK_R). Moreover, the core of our analysis was performed on the comparison between the Cd-treated samples of the two tissues indicated as Cd_R vs Cd L*. In this comparison, the Cd-treated root vs Cd-treated leaf DEGs (Cd_R vs Cd L*) were obtained by subtracting the DEGs belonging to the CK_R vs CK_L comparison (tissue specific DEGs in control conditions) to the Cd_R vs Cd_L. An adjusted p-value cutoff of 0.05 and a log 2 fold change (Log 2 FC) threshold of 1 was adopted to filter the significantly up-and down-regulated genes. A correlation analysis was performed in order to demonstrate experiment repeatability and to reveal differences in gene expression among samples. Principal Component Analysis 3D plot and a heatmap were obtained by using R language, considering the read counts of each sample, including biological replicates, as input data.

Real-time validation of selected DEG candidates using qRT-PCR
Total RNA (2.5 µg) purified from leaves and roots as described above, was reverse transcribed using Super-Script ™ Vilo ™ cDNA synthesis kit by ThermoFischer Scientific, according to manufacturer's instructions. Realtime qRT-PCR was carried out for a total of 8 DEGs with PowerUp SYBR Green Master mix by ThermoFischer Scientific in the Bio-Rad iQ5 Thermal Cycler detection system. All the genes have been normalized with cyclindependent kinase C-2 (CDKC-2, XM_004962139), being a suitable housekeeping gene [49]. All reactions were performed in duplicate and fold change measurements calculated with the 2 −∆∆CT method.

KEGG classification of cadmium and salt stress common genes
In order to elucidate the core response of A. donax to both salt [31,32] and cadmium stress, the significant DEGs (padj ≤ 0.05) belonging to both S3_vs_CK (specifically deregulated DEGs under severe salt stress, 256.67 mM NaCl), and S4_vs_CK (specifically deregulated DEGs under extreme salt stress, 419.23 mM NaCl) comparisons [31,32] were filtered considering a threshold of ± 2.00 log2FoldChange and merged, thus resulting in a "list of salt DEGs". The same procedure was adopted to retrieve the significant DEGs (padj ≤ 0.05) belonging to Cd_R vs Cd_L*, Cd_L vs Cd_CK and Cd_R vs Cd_CK comparisons, resulting in a "list of cadmium DEGs". A KO ID (KEGG ORTHOLOGY Database, https:// www. genome. jp/ kegg/ ko. html) was assigned to each DEG, then, the two lists were compared by using the KO ID as common annotation code. The aforementioned comparison generated a merged DEG list where only DEGs whose KO ID was in both the "list of salt DEGs" and "list of cadmium DEGs" were included. Finally, genes were grouped by concordant type of regulation (up-regulated or down-regulated in both stresses).

Effect of cadmium upon A. donax morpho-biometric and physiological parameters
As described in the Material and Methods section, both morpho-biometric and physiological parameters of A. donax G10 ecotype were evaluated at sampling date after being subjected to 4 ppm (4 mg/kg) of cadmium nitrate. A picture of giant reed G10 ecotype at sampling time is shown in Supplementary file 1: Figure S1. Considering the average values of the three biological replicates, we observed that both the main stem height and biomass dry weight per pot were reduced in those samples subjected to Cd treatment (Supplementary file 2: Figure S2A and S2B). Moreover, the net photosynthesis efficiency was also significantly reduced in treated samples compared to the untreated samples, reaching values of 8.64 µmol CO 2 m −2 s −1 and 14.91 µmol CO 2 m −2 s −1 , respectively (Supplementary file 2: Figure S2C). The alteration of the aforementioned parameters indicated the effectiveness of the cadmium dose to induce response in G10 giant reed ecotype.

Transcript assembly and annotation
In this study, the RNASeq approach was used to comprehensively identify the transcriptional response of the A. donax G10 ecotype in leaves and roots. A flow chart of the giant reed leaf and root transcriptome sequencing and de novo assembly process is shown in Supplementary File 3: Figure S3. RNA integrity was checked prior sequencing and the average RNA integrity number (RIN) was 7.05. This shows that all samples were of good quality for further processing and sequencing (Table 1). After sequencing, we filtered the raw reads to remove the adapter-based or poor quality reads, obtaining a total of clean reads equal to 416 million (Table 1) [29][30][31][32]. To assess assembly consistency, filtered unique reads were mapped to the reconstructed transcriptome and average read mapping rate using bowtie2 alignment software was equal to 69.0% (Table 1). The completeness of the assembled transcriptome was evaluated by comparing it to the set of Embryophyta genes using the BUSCO quality assessment tool coupled with the OrthoDB (9.0 version) database of orthologs [40]. Among the searched 1440 BUSCO groups, 68.54% (987 BUSCOs) were complete (912 singlecopy orthologs and 75 duplicated), 20.0% (288 BUSCOs) were fragment and 11.5% (165 BUSCOs) were missing. These results are comparable to those obtained previously in A. donax (70.07%, 13,19% and 16.74%, respectively) [29]. In addition, both transcript and unigene length distributions were reported (Supplementary file 4: Figure S4). Consequently, these results indicated that the sequencing quality was reliable to perform downstream analysis.  Table 2).

Identification of differentially expressed genes (DEGs)
The characterization of root and leaf A. donax transcriptome was accomplished by the identification of those unigenes whose expression level changed upon cadmium treatment. Based on the experimental design, a total of 162 genes showed differential expression in response to Cd treatment, 107 of them differentially expressed in root (Cd_R vs CK_R) and the remaining 55 differentially expressed in leaf (Cd_L vs CK_L). Among the 107 DEGs found in root, 36 were up-regulated and 71 downregulated. In leaf, 19 genes resulted up-regulated and 36 down-regulated (Fig. 1A). No DEGs were found in common between root and leaf. The higher number of DEGs in root with respect of the leaf suggested that root, representing the interface between soil and plant, is subjected to a wider reprogramming of the gene expression than leaves (Supplementary file 5; Table S1). DEGs identified in biological replicates clustered together in both organs, indicating good reproducibility of treatment (Supplementary file 6; Figure S5). Moreover, samples either belonging to control or treated samples clustered very close, thus probably explaining the low number of DEGs retrieved when comparing stressed and control samples (Supplementary file 6; Figure S5, S6). On the contrary, a total of 5303 genes were counted as differentially expressed when the Cd_R vs Cd_L* comparison was analyzed (Fig. 1B). This result was obtained by subtracting the DEGs found in the comparison CK_R vs CK_L (all genes that are  Figure S7). The results show high congruence between RNA-Seq and qRT-PCR (coefficient of determination R 2 = 0.91), which accounts for the high reliability of RNA-Seq quantification of gene expression.

Identification of transcription factor families involved in plant response to cadmium
Transcription factors (TFs) have been identified as candidate targets for ameliorating plant tolerance in case of cadmium treatment. Therefore, DEGs encoding for TFs were retrieved from each comparison (Cd_R vs CK_R, Cd_L vs CK_L and Cd_R vs Cd_L*), grouped according to the belonging family and sorted for their corresponding abundance (Fig. 4). As result of the analysis, the five more abundant transcription factor families resulted "Ethylene responsive transcription factor" (ERF) (29 DEGs), "bzip" (27 DEGs), "WRKY" (23 DEGs), "bHLH" (15 DEGs) and "MYB" (13 DEGs) (Fig. 4). When different comparisons were considered, two bHLH (both up-regulated) and two ERF (1 up-and 1 down-regulated) TFs were found in the Cd_R vs CK_R dataset and one ERF, two bZIP and one WRKY (all down-regulated) were found in the Cd_L vs CK_L dataset. On the contrary, a higher number of DEGs were found in the Cd_R vs Cd_L* dataset: 26 DEGs belonging to ERF family (21 up-and 5 down-regulated), 25 to bZIP family (14 up-and 11-down regulated), 22 to WRKY family (17 up-and 5 down-regulated), 13 to bHLH family (10 up-and 3 down-regulated) and 13 to MYB family (7 up-and 6 down-regulated).

Main Processes affected by cadmium treatment
To have a comprehensive view of the metabolic changes occurring in A. donax L. under Cd treatment, all the significant DEGs of the Cd_R vs Cd_L* comparison were mapped to the MapMan 3.5.1R2 pathways. The analysis indicated that the five most enriched metabolic pathways are "Protein biosynthesis", "Phytohormone action", "Nutrient uptake", "Cell wall organisation" and "Polyamine metabolism" ( Table 3). The DEGs belonging to these pathways resulted mostly up-regulated in treated root compared to the treated leaf (Cd_R vs Cd_L*), except for genes of "Polyamine metabolism" that resulted all down-regulated in root compared with the leaf. In addition, considering their crucial role in heavy metal stress response, the coding sequence of each transcripts belonging to "Oxidative stress", "Heavy metal ATPase 3 (HMA3) transporters" and "NRAMP transporters" categories along with the clusters resulted from the MapMan analysis (Table 3) have been aligned against the Arabidopsis thaliana genome (TAIR10) and the score of these alignments was reported (identity score and e-value) thus providing valuable indications of the cluster similarity with the reported genes (Table 4). Congruously, Table 4 reports clusters whose e-value was < 0.05 and showing a threshold of ± 2.000 log2fold change. The complete list of DEGs belonging to these pathways is reported in Supplementary file 9 (Table S2 A Regarding the "Phytohormone" category, nine-cisepoxycarotenoid dioxygenase 9 catalyzing the first step of abscisic-acid (ABA) biosynthesis from carotenoids is down regulated in cadmium treated giant reed roots (Cd_R vs Cd_L*) (Table 4). Moreover, abscisic acid 8'-hydroxylase 1 -related involved in the first step of the oxidative degradation of ( +)-ABA was up-regulated in the Cd_R vs Cd_L* comparison, suggesting that cadmium treatment suppresses ABA synthesis and promotes its degradation in giant reed root thus quenching the signal cascade induced by this hormone. Indole-3-pyruvate monooxygenase YUCCA2-related encoding an enzyme involved in auxin biosynthesis and converting the indole-3-pyruvic acid (IPA) to indole-3-acetic acid (IAA) resulted up regulated in the cadmium treated roots with respect to treated leaf (Cd_R vs Cd_L*). Similarly, the gene encoding 1-aminocyclopropane-1-carboxylate oxidase (ACO) encoding an enzyme synthesizing the plant hormone ethylene was up-regulated in giant reed roots treated with cadmium. Gibberellin 20 oxidase 2, encoding a key oxidase enzyme in the biosynthesis of gibberellins, considered potent signaling molecules to enhance growth and developmental processes under abiotic stress conditions [52], was also found up-regulated in the Cd_R vs Cd_L* comparison. Signaling peptides, also known as peptide hormones, along with classical phytohormones, are involved in plant intracellular signaling. C-terminally encoded peptide (CEP), two genes encoding ROTUNDIFOLIA like proteins, and a transcript encoding CLAVATA3/ESR-RELATED 25 were among the sharply up-regulated genes in the cadmium treated root (Cd_R vs Cd_L*).
As concerns the "Transcription factors" category, several transcripts encoding Ethylene-Responsive Transcription (ERT) factors which have been reported to be involved in the ethylene signaling transduction pathway in plant abiotic stress response [53] were found sharply up regulated in the Cd_R vs Cd_L*. Therefore, the results indicate that ethylene-triggered signal cascade plays a crucial role in giant reed root under cadmium treatment. Interestingly, the WRKY gene family is also over-represented in the "Transcription factor" category (Table 4) and in particular WRKY9, that is involved in increased root suberin deposition [54], WRKY33 that controls the apoplastic barrier formation in roots [55] and WRKY50 positively regulating resistance against necrotrophic pathogens [56] were found up-regulated in cadmium treated giant reed roots (Cd_R vs Cd_L*). Conversely, WRKY46 that seems to be involved in hypersensitivity to drought and salt stress [57] was down regulated ( Table 4).
The analysis of the "Cell wall" category reveals that transcripts encoding alpha expansin 11, that causes loosening and extension of plant cell walls by disrupting non-covalent bonding between cellulose microfibrils and matrix glucans [58], were up-regulated in the Cd_R vs Cd_L* comparison. The genes encoding cinnamoyl Coa reductase, laccase and membrane-associated progesterone binding protein 3 involved in lignin biosynthesis were strongly induced in cadmium treated roots. Moreover, pectin methylesterase 61 and invertase/pectin methylesterase inhibitor gene have been found upregulated in Cd-treated root thus confirming the overall results related to "Cell wall" category that this cellular district is target of deep remodeling under cadmium stress [59]. Regarding the "Nutrient uptake" category, transcripts related with nitrogen metabolism, such as nitrate reductase and high affinity nitrate transporter 2.5 were upregulated by cadmium treatment in giant reed roots, being they involved in nitrogen assimilation and uptake, respectively. Furthermore, phosphate transporter encoding genes involved in phosphate uptake at plasma membrane level were also upregulated whereas vacuolar iron transporter required for iron sequestration into vacuoles was downregulated in cadmium treated roots. Overall, these results support that cadmium treatment induced an imbalance of the main nutrient levels which can cause disturbances in plant growth. Polyamines (PAs) are ubiquitous low molecular weight aliphatic cations that are present in all organisms; the major PAs in plants are putrescine, spermidine, and spermine, and to a lesser extent, cadaverine. Cadaverine, a structurally different diamine, has an independent biosynthetic pathway as it is synthesized from lysine by lysine decarboxylase (LDC) [60]. Agmatine deiminase catalyzes the hydrolysis of agmatine into N-carbamoylputrescine in the arginine decarboxylase (ADC) pathway of putrescine biosynthesis. In the "Polyamine metabolism" category, lysine decarboxylase gene was upregulated whereas agmatine   deiminase was downregulated in cadmium treated roots, suggesting that under heavy metal stress, the pathway leading to cadaverine is preferentially undertaken in giant reed root. The amount of cadmium inside the root cells depends on the balance between metal uptake, normally determined by general metal transporters belonging to the NRAMP family, and metal extrusion. The heavy metal ATPase 3 (HMA3) transporters translocate cadmium both outside the cell and into the vacuole [61,62]. The analysis of "Metal transporter" group revealed that transcript encoding NRAMP1 transporter was upregulated whilst HMA3 transporter encoding gene were downregulated in the cadmium treated root thus indicating that the mechanisms to avoid cadmium uptake and stay are not implemented in our conditions.
As concerns the ROS scavenging regulatory mechanisms, the analysis of Cd_R vs Cd_L* revealed that superoxide dismutase (SOD), ascorbate peroxidase (APX), catalase and NADPH oxidase 1 encoding genes are all up-regulated suggesting that H 2 O 2 could be the main ROS the A. donax roots have to cope with under cadmium treatment. Moreover, many transcripts encoding glutathione transferases (GSTs) are also up-regulated (Table 4) concordantly with their role in abiotic stress relief [63,64]. Moreover, transcripts encoding glutaredoxin-C5, small redox enzyme of approximately one hundred amino-acid residues involved in the reduction of specific disulfides such as glutathionylated proteins [65] were found strongly induced in cadmium treated roots. Finally, the "Protein biosynthesis" category encompasses a plethora of genes involved in ribosome constitution (ribosomal proteins) and DNA-directed RNA polymerase II subunit that resulted up regulated indicating that this pathway is sharply activated by cadmium treatment in giant reed roots (Table S2C).

Analysis of the A. donax response to salt and cadmium treatments
The availability of transcriptomic data enabled the identification of genes that are implicated in A. donax response to both salt [31,32] and cadmium stress. In details, a DEG list containing deregulated genes in both stress conditions was obtained by using the KEGG ORTHOL-OGY database (see Materials and Methods). The KEGG ORTHOLOGY database represents a reference resource for gene and protein annotation and allowed to overcome the species-specific annotation of other databases potentially leading to bias in the analysis [66]. Figure 5 reports the classification of common DEGs and their numerical distribution. Based on this analysis 162 DEGs were successfully classified and most of them (156) were up-regulated in both stress conditions whereas only 6 DEGs were down-regulated. In particular, the most abundant KO category was (ko01000) Enzymes (51 up-regulated and 4 down regulated DEGs), followed by (ko03011) Ribosome (46 up-regulated DEGs), (ko04147) Exosome (17 up-regulated DEGs), (ko04131) Membrane trafficking (15 up-regulated DEGs) and (ko02000) Transporters (13 upregulated and 1 down-regulated DEGs) (Fig. 5). Among the DEGs included in the "Transcription factor" category AP2-like factor and EREBP-like factor were found upregulated, along with amino-cyclopropane carboxylate oxidase (ACO) confirming the crucial role of ethylene in giant reed abiotic stress response. Moreover, the "Ribosome" category (46 up-regulated DEGs) contains an elevated number of genes encoding ribosomal proteins indicating that "Protein biosynthesis" is an extremely involved pathway in both conditions. The complete list containing the KO description of each DEG involved in both salt and cadmium stress conditions is in Supplementary file 10: Table S3.

Discussion
Anthropogenic activities such as metal industries, mining application of pesticides and fertilizers led to dangerous Cd accumulation in arable soil worldwide [61]. The high solubility of Cd in the soil associated with the plant capability to absorb it represents a threat to humans as final consumers of putatively contaminated fruit and vegetables. The strategy to allocate cadmium contaminated soil to bioenergy crop might turn out to be successful as it potentially solves the increasing demand of sustanible energy sources which can be achieved without impairing the quote of agricultural lands. Arundo donax L. is considered tolerant to several abiotic stress and recently, due to its ability to accumulate high concentration of HMs [24,25] has been proposed as a potential phytoremediation species. In order to elucidate the response of giant reed to cadmium, in this work we sequenced and de novo assembled the A. donax L. leaf and root transcriptome after a prolonged soil treatment with 4 mg/Kg cadmium nitrate. Cadmium concentrations in uncontaminated soil usually is 0.5 mg/Kg and cadmium contents above 3 mg/kg are generally thought to indicate contaminated soil [37]. Based on the measurement of trace cadmium concentration under non-phytotoxic cadmium treatment, Bonanno (2012) [67] showed that giant reed roots act as main accumulation centre (root to shoot ratio is 6.6:1), whereas stem functions as transition centre given the involvement in ion translocation from root to leaves. Taken into account these observations which preliminarily indicated that root is the frontline in encountering potentially toxic cadmium levels, we focused our analysis on a Cd_R vs Cd_L* comparison originated by subtracting the DEGs belonging to the CK_R vs CK_L comparison (root specific DEGs under control conditions) to the Cd_R vs Cd_L comparison (root specific DEGs under cadmium treatment) thus obtaining a pool of DEGs specifically up or down-regulated (compared to leaf DEGs) in treated roots. Based on differentially gene expression data, the up regulation of the NRAMP1 transporter, a crucial group of transmembrane protein involved in the uptake of a broad range of ions (Mn, Zn, Cu, Fe, Ni, Co and Cd), suggests that giant reed roots reinforce their ability to adsorb ions in the presence of cadmium which interferes with the essential mineral supply [61]. Moreover, the down regulation of HMA3 transporter encoding gene, that regulates cadmium efflux out the citosol towards the vacuole or outside the cell [51], indicates that the mechanisms to avoid cadmium residence in root cells are not implemented in our conditions. Once in the citosol, Cd has to be detoxified in order to avoid the onset of cellular damage. It has been shown that the major mechanism of cadmium detoxification is based on phytochelatins, a family of cysteine-rich oligopeptides synthesized from glutathione (GSH) by the enzyme phytochelatin synthase (PCS), that chelates the HMs to form complexes readily stored in the vacuoles [68]. Recently, three giant reed PCS genes have been characterized both in the species of provenance and in transgenic model organisms [34]. These genes, namely AdPCS1-3, are expressed in all organs (root, rhizome, node, internode, leaf sheath and blade) in normal condition. In response to 500 μM CdSO 4 treatment, all the PCS genes showed an increase of expression in roots whereas any change in leaf gene expression was observed [34]. Although these results clearly indicate that A. donax genome contains at least three PCS genes constitutively expressed in all plant organs and whose expression is induced in roots under cadmium treatment, PCSs were found among the DEGs (data not shown) suggesting that the detoxification process through metal chelation was not triggered in our conditions. After the stress perception, phytohormones are the chemical messengers that play a pivotal role in the induction and regulation of diverse signal transduction pathways in response to cadmium stress [69]. Higher endogenous abscisic acid (ABA) concentration are reported to mitigate the damaging effects of Cd stress by promoting root-to-shoot Cd translocation through the apoplast and more Cd accumulation in the shoots [70]. Our results suggested that abscisic acid biosynthesis is inhibited whereas the main gene involved in ABA degradation is up-regulated in giant reed Cd-treated roots (Table 4), probably causing Cd to be detained in the roots. Auxin is recognized as a crucial phytohormone in regulating every aspect of plant root development during normal and stress conditions [71]. In Arabidopsis thaliana, significantly enhanced expression of auxin biosynthesis gene YUCCA6 was observed upon Cd treatment [72]. Moreover, reduced transcript levels of auxin efflux carrier pin-formed 1 (PIN1) were detected under Cd exposure in the post-embryonic roots decreasing auxin transport in the root apex thus altering auxin homeostasis [72]. In our study, both auxin biosynthesis and transport transcript (YUCCA2 and PIN) were found among the up-regulated genes in Cd-treated giant reed roots suggesting that roots attempt not to modify auxin homeostasis preventing the alteration of root architecture as observed in tobacco [73]. Furthermore, our results strongly indicate that ethylene biosynthesis and the downstream signaling cascade are up regulated in Cd-treated roots (Table 4). Similarly, the enhanced transcript levels of 1-aminocyclopropane-1-carboxylate (ACO) and ethylene signaling genes (EIN) were observed under As applicationin the rice roots [74], and transcript level accumulation of ethylene receptor 2 (ETR2) gene was also observed under Cd application in Arabidopsis thaliana generally suggesting the occurrence of Cd-mediated induction of ethylene production and signaling [75].
Several studies indicated that ethylene modulates the ROS machinery by regulating the activities of both ROS producing and scavenging enzymes, this being considered the crucial reaction following HM stress in plants [76]. Ethylene stimulates the production of ROS by activating nicotinamide adenine dinucleotide phosphate (NADPH) oxidases under HM toxicity [77] and also modulates antioxidant defense systems in plants exposed to HM stress. NADPH oxidase 1 transcript was among the up regulated genes in Cd-treated roots and the activation of catalase, SOD, APX and GST expression (Table 4) is probaly an attempt to overcome the Cd-induced oxidative stress. This result assumes an important significance since in several species it has been reported that the activation of ROS scavenging machinery is achieved in ethylene insensitive mutants [69] thus suggesting that this type of response could be specific of A. donax roots. However, we cannot exclude that the induction of YUCCA2 involved in auxin biosynthesis might be responsible of the regulation of reactive oxygen species (ROS) homeostasis as showed in transgenic potato over-expressing YUCCA6 gene [78]. Oxidative stress can induce protein S-glutathionylation modulating protein function and to prevent irreversible oxidation of protein thiols [79]. High levels of oxidized glutathione (GSSG) can be sufficient to trigger protein S-glutathionylation by a thiol-disulphide exchange reaction between a protein thiol and GSSG. Interestingly, transcript encoding glutaredoxin-C5 is sharply up regulated in giant reed roots indicating that an active reactivation of protein function via the reduction of the cysteine sulphydrilic groups might occurr.
Within the "Phytohormones" category, two genes encoding ROTUNDIFOLIA like proteins, which are small polypeptides acting as a regulatory molecules coordinating cellular responses involved in different aspects of cell differentiation, growth and development, were found sharply up regulated in the Cd_R vs Cd_L* comparison. Their role is uncertain but probably they act by restricting polar cell proliferation in lateral organs and coordinating socket cell recruitment and differentiation at trichome sites [80]. A transcript encoding CLAVATA3/ ESR-RELATED 25 is among the sharply up-regulated genes in the cadmium treated root (Table 4), CLAV-ATA3/ESR represents a group of plant proteins acting as extracellular signal peptides involved in cell-to-cell communication in concert with different receptors in a range of processes during plant development [81].
A transcript encoding C-terminally encoded peptide (CEP) was strongly up-regulated in cadmium treated roots. CEP is a 15-amino acid post-translationally peptide which plays a pivotal role in lateral root formation and nodulation and its overexpression in Medicago results in the inhibition of lateral root formation and enhancement in nodule formation. Besides the role in root architecture, CEP genes also play a role in nitrate assimilation under N starvation conditions and provoke the up regulation of the nitrate transporter gene NRT2.1 in roots specifically when nitrate is present in the rhizosphere [82]. Nitrate transporter are involved in the constitutive high-affinity transport system (cHATS) under low nitrate conditions. The principal role of this cHATS is to enable roots previously deprived of nitrate to absorb this ion and initiate induction of nitrate-inducible genes [83]. Cd has been reported to inhibit NO 3 − uptake in several plant species, under normal and high nitrate supply [84]. Therefore, it is possible that nitrate uptake is inhibited under Cd stress in giant reed roots mimicking low-nitrate conditions. Interestingly, the analysis of "Nutrient uptake" category revealed that nitrate transporter 2.5 genes are strongly up-regulated (Table 4) suggesting that mechanisms to enhance nitrogen supply are implemented in A. donax cadmium treated roots. Moreover, nitrate reductase encoding gene, involved in nitrate assimilation by reducing nitrate to nitrite, is up-regulated in cadmium treated samples probably to provide for more enzyme molecules that are supposed to be inhibited by cadmium [85]. Under cadmium stress a decrease of mineral nutrient concentrations in plant leaf, such as Fe and P, has been observed and represents the key reason for the restraint of leaf photosynthesis [61]. In our work, transcript encoding the plasma membrane phosphate transporter are up-regulated in Cd-treated roots probably to face the unavailability of soluble phosphate sequestered in the soil as cadmium phosphate [86].
The plant cell has a variety of mechanisms tools to avoid Cd stress. Mainly, cell wall remodeling might prevent Cd from entering and damaging the protoplast [87]. At primary cell wall, pectin, which contains most of the negative charges, can immobilize Cd very effectively. Furthermore, secondary cell wall lignification can serve to create a barrier to prevent cadmium entry. The utilization of different antibodies detecting methylated and demethylated forms of pectin in cadmium stressed Linum usitatissimum hypocotyls has led to the identification of low-methylated pectin, which is particularly able to bind Cd ions due to the presence of free carboxyl groups in the outer side of the primary cell wall. The increased low-methylated pectin form occurs along with the colocalized upregulation of pectin methylesterase (PME) activity. Instead, a higher amount of methylated pectin was detected within the inner side of the primary cell wall which was hypothesized to have an impermeable function aimed to keep the cytosolic Cd away from the cell [88]. In this perspective, the up regulation of both PME and PME inhibitor encoding genes might serve to finely regulate Cd sequestering at primary cell wall in giant cane root. At secondary cell wall level, lignification makes the cell wall less penetrable thus creating an effective barrier against the entry of Cd [89]. The induction of the lignification process appeared as a key process useful to discriminate Cd-sensitive and -tolerant plants [90,91]. The discovery of the upregulation of several WRKY transcription factors involved in cell wall lignification as well as the induction of cinnamoyl CoA reductase, laccase and membrane-associated progesterone binding protein 3, all of them involved in lignin biosynthesis (Table 4), clearly indicate that giant cane roots respond to cadmium treatment by avoiding its entrance into the cell.

Conclusions
The global analysis of our findings suggests that prolonged cadmium exposure stimulated a clear response at both morpho-physiological and transcriptomic levels. Hence, cadmium treated plants showed significantly reduced main stem height, biomass dry weight and the net photosynthesis efficiency. The quality of transcriptome sequencing and assembly was elevated and led to the identification of crucial metabolic pathways and to decipher the A. donax response to cadmium stress. Three main factors have to be taken in strong consideration in this concluding remarks: a) the used cadmium concentration (4 mg/Kg), slightly higher than allowed; b) the induction of lignification process clearly suggested by transcriptome analysis; c) the lack of phytochelatin transcripts among the DEGs. In our opinion, all these issues indicate that the induced stress condition can be sensed as "mild stress". The low number of DEGs within the CK_R vs Cd_L and CK_L vs_Cd_L comparisons seems to be in line with this hypothesis. The undertaken strategy was to analyse the Cd_R vs Cd_L* comparison and it led us to focus on the main patterns involved in the Cdtreated giant cane root, it being the interface between plant and soil. Our results suggest that ethylene biosynthesis and signaling are strongly activated. In this respect, the identification of several genes differently regulated in both salt and cadmium conditions, such as genes involved in ethylene biosynthesis and signal transduction, outlines those metabolic pathways and biochemical reactions as useful markers of abiotic stress in giant reed. Finally, the finding of DEGs encoding several small peptides functioning as messenger molecules between root and shoot in order to communicate the stressful status to the upper part of the plants (CEP, ROTUNDIFOLIA), the induction of the ROS scavenging machinery, and, above all, the remodelling of plant cell wall confirm the tolerance of giant cane towards cadmium stress and strongly support its cultivation in cadmium contaminated soils in a perspective to save agricultural soil for food and feed crops.