Differential levels of gene expression and molecular mechanisms between red maple (Acer rubrum) genotypes resistant and susceptible to nickel toxicity revealed by transcriptome analysis

Abstract Knowledge of regulation of genes associated with metal resistance in higher plants is very limited. Many plant species have developed different genetic mechanisms and metabolic pathways to cope with metal toxicity. The main objectives of this study were to 1) assess gene expression dynamics of A. rubrum in response to nickel (Ni) stress and 2) describe gene function based on ontology. Certified A. rubrum genotypes were treated with 1,600 mg of Ni per 1 Kg of soil corresponding to a soil total nickel content in a metal‐contaminated region in Ontario, Canada. Nickel resistant and susceptible genotypes were selected and used for transcriptome analysis. Overall, 223,610,443 bases were generated. Trinity reads were assembled to trinity transcripts. The transcripts were mapped to protein sequences and after quality controls and appropriate trimmings, 66,783 annotated transcripts were selected as expressed among the libraries. The study reveals that nickel treatment at a high dose of 1,600 mg/kg triggers regulation of several genes. When nickel‐resistant genotypes were compared to water controls, 6,263 genes were upregulated and 3,142 were downregulated. These values were 3,308 and 2,176, respectively, when susceptible genotypes were compared to water control. The coping mechanism of A. rubrum to Ni toxicity was elucidated. Upregulation of genes associated with transport in cytosol was prevalent in resistant genotypes compared to controls while upregulation of genes associated with translation in the ribosome was higher in susceptible genotypes when compared to water. The analysis revealed no major gene associated with Ni resistance in A. rubrum. Overall, the results of this study suggest that the genetic mechanism controlling the resistance of this species to nickel is controlled by genes with limited expression. The subtle differences between resistant and susceptible genotypes in gene regulation were detected using water‐treated genotypes as references.


| INTRODUC TI ON
Plants deal with toxic levels of metals such as nickel (Ni) using different mechanisms. Ni is a transition metal that is found in natural soils at low concentrations. An excess of Ni ions in the soil results in the disorder of cell membrane functions, inhibition of cell division in the root system, and a decrease of nontolerant plant growth (Yadav, 2010). It can also increase the concentration of hydroxyl radicals, superoxide anions, nitric oxide, and hydrogen peroxide (Bhalerao, Sharma, & Poojari, 2015;Boominathan & Doran, 2002;Rao & Sresty, 2000;Stohs et al., 2001). Metal ions that accumulate in plant systems usually disturb cellular ion homeostasis and can generate (directly or indirectly) reactive oxygen species (ROS) which may cause oxidative stress. The toxicity of ROS can lead to the destruction of DNA structure and enhance the oxidation of lipids and proteins.
Nickel has a complex chemistry which complicates the decryption of its toxicity mechanisms in plants (Bhalerao et al., 2015). It does not directly induce the production of ROS as it is not a redoxactive metal. Its role in ROS production is an indirect one by inhibiting the function of several antioxidant enzymes which include ascorbate peroxidase (APX), catalase (CAT), glutathione peroxidase (GSH-Px), glutathione reductase (GR), guaiacol peroxidase (GOPX), peroxidase (POD), and superoxide dismutase (SOD) (Bhalerao et al., 2015;Freeman et al., 2004;Gomes-Junior et al., 2006;Pandey & Sharma, 2002).
Nickel translocation and toxicity in hardwood species has been recently a focus of a few reports. It has been demonstrated that white birch (Betula papyrifera), trembling aspen (Populus tremuloides), and red oak (Quercus rubra) accumulate nickel in leaves. Therefore, they are classified as nickel accumulators (Mehes-Smith & Nkongolo, 2015;Theriault, Michael, & Nkongolo, 2016a;Theriault, Nkongolo, & Michael, 2014;Theriault et al., 2013;Tran et al., 2014). Red maple (Acer rubrum) on the other hand does not accumulate nickel in its tissues (the amount of bioavailable nickel in the soil is higher than the total nickel in roots). The translocation of nickel from roots to aerial parts is also very small. This species can be therefore classified as a nickel avoider (Kalubi, Mehes-Smith, & Omri, 2016;Kalubi et al., 2015). A. sacharinium, a close relative of A. rubrum, stores Ni in its roots with limited translocation to other plant parts. It is classified as a Ni excluder . We anticipated that genetic mechanisms and metabolic pathways involved in coping with nickel toxicity in these hardwood species, with distinct strategies in dealing with this element in their system, would be different.
Most resistant plants that accumulate metals in their tissues have developed a detoxification mechanism that involves phytochelatin (PC). These tri-peptides are synthesized from reduced glatathione (GSH). GSH conjugates with heavy metal molecules through glutathione S-transferase during detoxification processes. Many studies have demonstrated that genes controlling glyoxalases, phytochelatin synthase, glutathione reductase, serine acetyltransferase, ATP sulfuylase, cystathionine synthase, glutathione synthetase, and γglutamycylcysteine synthetase are candidates for providing metal tolerance by regulating GSH and PCs levels in accumulator plants (Yadav, 2010). Our understanding of genetic mechanisms involved in metal avoidance in plant species such as A. rubrum is vague.
The main objectives of this study were to (1) assess gene expression dynamics of A. rubrum in response to nickel stress and (2) describe gene function based on ontology. The study provides the first description of molecular and biological processes involved in Ni resistance in A. rubrum.

| Nickel treatment
Acer rubrum seeds were provided by the National Tree Seed Centre, Canadian Forest Services (New Brunswick, Canada). These certified seeds (accession # 2001 1031.0) were collected from Larry Brook, New Brunswick (Canada). Assessment of nickel toxicity is described in Theriault and Nkongolo (2016) and Theriault et al. (2016a). The experimental layout was a completely randomized design with one Ni treatment and two types of control. To assess the genetic resistance of A. rubrum genotypes to Ni, 45 six-month seedlings were treated with 1,600 g of Ni per 1 kg of dry soil using nickel nitrate [Ni(NO 3 ) 2 ] salts and grown in a growth chamber. This concentration that was used in previous studies corresponds to the level of total nickel in contaminated sites in the mining region of the City of Greater Sudbury. Details of seed germination and seedlings treatment with nickel nitrate are presented in Theriault and Nkongolo (2016) and Theriault et al. (2016a). Genotypes treated with water only were used as the main control. Potassium nitrate (KNO 3 ) treatment was used to control for the nitrate effects. The treatment and the controls were replicated 15 times. Damage rating was recorded every two days based on a scale of 1 to 9, 1 = no visible toxicity symptoms and 9 = dead plants. Individual plants with a score of 1 to 3 were considered nickel resistant, 4 to 6, moderately resistant, and 7 to 9 susceptible (Theriault & Nkongolo, 2016). Genotypes resistant and susceptible to a soil nickel concentration of 1,600 mg/kg are analyzed in detail. For transcriptome analysis, three Ni resistant and three Ni susceptible genotypes were selected along with three genotypes from each of the controls (water and nitrate controls).

| De novo transcriptome assembly
Methods for extraction, RNA-seq libraries, next-generation sequencing, and de Novo transcriptome assembly are detailed in Theriault, Michael, and Nkongolo (2016b). The libraries were quantified using Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA) and the sequencing was performed on the Illumina HiSeq 2000 sequencing system (Illumina Inc.) at Seq Matic (Fremont California, USA).
The RNA-seq data from all the samples including six nickel-treated (three resistant and three susceptible), three water-treated (control), and three nitrate-treated were used as input for the Trinity program (http:trinityrnaseq.githb.io) to assemble the transcripts. The raw reads were mapped to Trinity assembled transcripts using bowtie (http://bowtie-bio.sourceforge.net/index.shtml), and RSEM (http:// deweylab.biostat.wisc.edu/rsem) was used to quantify transcript and expression levels. Additional QC at transcript level was performed, including number of unique transcripts detected, percentage of reads belonging to the top transcripts expressed, normalization for RNA composition, and grouping, and correlation between samples.
If a transcript had a count per million (CPM) value ≥1 in at least two of the samples, we considered it expressed in the experiment and included it for downstream QC analysis. Gene expression was calculated and expressed as Reads Per Kilobase per Million reads mapped (RPKM) (Mortazavi et al., 2008). The count per million (CPM) cutoff was 0.63 based on the average read count of all samples (15.9 million). This CPM cutoff roughly equaled to 10 raw reads in this experiment. A gene with a CPM value >0.63 in at least two samples from the experiment was included for downstream analysis.
Multidimensional plots were created to view sample relationships. This was performed using R Limma package. We also used the made4 (multivariate analysis of microarrays data using ADE4) program to cluster samples and drew heatmaps based on genes that had variable expression across samples. These variable genes were chosen based on a standard deviation (SD) of expression values larger than 30% of the mean expression value. Genes with a mean logCPM <1 were removed.
For differentially expressed genes, the normalized data were transformed to log2CPM values using the voom method from the R Limma package. A linear model was built for each comparison using the R Limma package, and statistics for differential expression analysis were calculated. The statistical values included log fold change (logFC), p-value, and false discovery rate (FDR). An FDR of 0.05 was used as the standard cutoff (two-fold change) to determine differentially expressed genes between treatments.
All transcripts were mapped to protein sequences in the UniProt database (http:// www.uniprot.org/) and the best match was used to annotate genes and assign gene ontology information. The annotated sequences were run through the GO-Slim function of the BLAST2GO program to provide a summary known gene functions (Conesa & Götz, 2008). The ontology categories included biological process, cellular components, and molecular function.

| Acer rubrum resistance to nickel
The genotypes screened in this study expressed a high level of resistance to Ni. In fact, most seedlings treated with nickel nitrate showed no damage (rating of 1 on the 1 to 9 scale) after two weeks of treatment with the exception of three genotypes that were susceptible (rating of 7). No plant damage from the nitrate control treatment was observed. Transcripts were assigned ontology and grouped by biological process, molecular functions, and cellular components. For biological processes, 58% of the 15,078 transcripts that were assigned ontology were identified under the following categories:
For molecular functions, 35% of transcripts code for proteins involved in nucleotide binding activities, 12.05% kinase activities, 11% DNA binding, and 10.2% transport activities ( Figure 2). For cellular component, 21.4% of the 13,676 transcripts that were assigned ontology were localized in the cytosol, 19.39% in the ribosome, 8.34% in the endoplasmic reticulum, 7.44% in the cytoskeleton, 7% in the plasma membrane, and 5.5% in the Golgi apparatus ( Figure 3).
Among the three principal ontologies, most of the expressed transcripts were classified into cellular component organization, transport, carbohydrate metabolic process, translation, catabolic process, and response to stress suggesting that these functional processes play a major role in A. rubrum gene activities.

| Differential gene expression
Hierarchical clustering can provide good indications of sample and gene relationships. The overall heatmap shows that clustering was good with all the controls (water and nitrate controls were similar to each other). Any differential gene expression was attributed to nickel treatments. Hence, water and nitrate control data were pooled. After normalization, a total of 66,783 transcripts were selected as differentially expressed. Surprisingly, there were no significant differences F I G U R E 2 Gene ontology of transcripts from red maple (Acer rubrum) control plants (water only). A total of 13,676 transcripts were grouped under molecular function using the BLAST2GO software. * TFA stands for Transcription factor activity

| Pairwise comparison of resistant and susceptible genotypes to water controls
The top upregulated and downregulated transcripts from the differential expression analysis were ranked based on logFC (Tables 1-4 and Tables S1-S4). The study reveals that the nickel treatment at a Pairwise comparisons for biological process, molecular functions, and cellular compartments between RG and water and SG and water are described in Figures 4-6. For biological process, the most differentially expressed transcripts coded for proteins associated with transport, cellular component organization, catabolic process, carbohydrate metabolic process, response to stress, translation, lipid metabolic process, and cellular modification process when RG were compared to water (Figure 4). There was twice the number of upregulated compared to downregulated transcripts that coded for transporters and proteins involved in translation activities in RG compared to water control. A small proportion of differentially expressed transcripts translated proteins associated with flower development, response to biotic stimuli, secondary metabolic process, cell growth, embryo development, and regulation of gene expression and epigenetics. For molecular function, over 30% of the transcripts that were upregulated or downregulated coded for proteins that were associated with nucleotide binding ( Figure 5). There was more downregulated (17%) than upregulated (9.5%) transcripts coding for proteins associated with kinase activity when RG were compared to controls. For cellular component, most differentially expressed transcripts coded for proteins found in the cytosol. The majority of these transcripts were upregulated as there were three times more upregulated than downregulated transcripts in that category. A different trend was observed in the cytoskeleton, plasma membrane, plastid, with more downregulated than upregulated. The transcripts coding for proteins found in the cell wall, vacuole, and thylakoid were all downregulated in RG when compared to water ( Figure 6). Hence, upregulation of genes associated with transport in cytosol is the main mechanism involved in RG in the presence of Ni.
When SG were compared to water, the top biological processes were translation, cellular component organization, transport, carbohydrate metabolic process, response to stress, catabolic process, and signal transduction (Figure 7). There were more upregulated transcripts than downregulated for translation and signal transduction while no difference between upregulation and downregulation was observed for genes associated with the other main biological processes. A downregulation was observed for transcripts coding for proteins associated with lipid metabolic processes and anatomical structure morphogenesis. As with RG, few downregulated transcripts coded for proteins associated with flower development, secondary metabolic process, embryo development, and regulation of gene expression and epigenetics.
For molecular function, the pattern of gene regulation was similar to that observed when RG were compared to water (Figure 8).
Transcripts coding for proteins associated with nucleotide bind-

| Genetic resistance to nickel
Because of their lack of mobility, plants are continuously exposed to abiotic stresses.

| Gene expression and ontology analyses
We found no difference in gene expression at high stringency (FDR >0.05) between Ni resistant and susceptible genotypes. The heatmap clustering also showed that the two types of genotype (RG and SG) form a distinct cluster that was separated from the control groups.
RG and SG samples intermingled within the same cluster. Some differences in gene expression were observed at low stringency based on a p-value analysis. But, p-value is more susceptible to including false-positive results than the False Discovery Rate. Hence, the results based on the p-values analysis were not considered in the analysis of gene expression. Significant differences between RG and SG were found when they were compared to water. The validation of DEG results by qPCR was not necessary because a series of quality controls during assembly and sequencing, and data processing was performed. In addition, data were analyzed using stringent statistical tests. Moreover, we used two types of references (water and nitrate controls) to filter the effect of nitrate on gene expression.
Furthermore, the amount of RNA recovered specially in susceptible genotypes was not enough to run other types of validation analysis such as qPCR.
In general, gene ontology (GO) revealed relevant information on possible mechanisms involved in Ni resistance in A. rubrum. GO defines gene functions and how these functions are related to each other. It describes molecular function (molecular-level activities performed by gene products), biological process (the larger processes, or 'biological programs), and cellular components (the locations relative to cellular structures in which a gene product performs a function) (Gene Ontology Consortium, 2017). In the present study, there was an upregulation of genes associated with transport in cytosol in RG compared to water control while upregulation of genes associated with translation in the ribosome was prevalent in susceptible genotypes when compared to the same control. In contrast, Theriault et al. (2016b) reported that the main mechanism involved in nickel resistance in B. papyrifera, a Ni accumulator, is a downregulation of genes associated with translation (in ribosome), binding, and transporter activities.
Moreover, in B. papyrifera, six candidate genes associated with nickel resistance were identified (Theriault et al., 2016b). They in- Overall, the complexity of the mechanism involved in nickel resistance in A. rubrum could be associated with the physiological process used by this species to cope with nickel contamination. The metal avoidance in plants such as A. rubrum is associated with morphological changes at the root system and likely involves auxins (Cai et al., 2011;Khare et al., 2017;Liu et al., 2011;Overvoorde, Fukaki, & Beeckman, 2010;Potters et al., 2007;Vitti et al., 2014). Plants also limit metal assimilation by the roots by secreting a number of substances such as organic acids, and substances in root extracellular matrix such as sugars, phenols, amino acids, and polysaccharides (Cai et al., 2011;Guo, Liang, & Zhu, 2009;Jutsz & Gnida, 2015). The mechanism used by A. rubrum to avoid nickel is unclear, but considering that the nickel treatment was conducted in a controlled environment using a sterilized sand / soil mix, it is likely that the avoidance mechanism is in situ.  (Meharg, 2005). Subcellular and intracellular localization are unlikely in A. rubrum considering the absence of Ni accumulation in its tissues. The lack of expression of glutathione associated genes that play a key role in redox stress processes suggests that this mechanism is not involved as well in A. rubrum response to nickel contamination (Hartley-Whitaker, Ainsworth, & Meharg, 2001;Hartley-Whitaker, Woods, & Meharg, 2002;Meharg, 2005). Our data suggest that uptake kinetics and metabolisms are the key processes involved in A. rubrum reaction to Ni toxicity. Mechanisms for protecting plants from oxidative stress appear to be constitutive in both resistant and sensitive genotypes.
The results of the present study suggest that resistance genes preexist in the genome of A. rubrum considering that the seeds used for this investigation are from an A. rubrum population that was not previously exposed to Ni. Other studies have reported Ni resistance in genotypes from metal-contaminated sites (Kirkey, Matthews, & Ryser, 2012;Watmough & Hutchinson, 1997). The fact that no major genes associated with metal transport or processing were identified in the analysis of the RG and SG A. rubrum transcriptome suggests that the genetic mechanism controlling the response of this species to nickel is controlled by several genes with limited expression.
The subtle differences between resistant and susceptible in gene regulation, if they exist, were difficult to detect in direct pairwise comparison.

| CON CLUS ION
The main goal of the present study was to determine the mechanism involved in A. rubrum response to nickel toxicity. Nickel treatment at a high dose of 1,600 mg/kg triggers regulation of several genes. The study revealed that unlike B. papyrifera, there were no significant differences in genes expression when RG and SG were compared based on the FDR test. However, distinctive differences between the two groups were found when they were compared to water controls. Upregulation of genes associated with transport in cytosol is the main mechanism associated with response to Ni in RG while upregulation of genes involved in translation in the ribosome was prevalent in susceptible genotypes in comparison with control. No major genes associated with nickel resistance were identified. Several highly expressed genes were found in the top 100 upregulated and downregulated based on heatmap profiles and logFC analysis. This suggests that the genetic mechanism controlling the resistance of A. rubrum to nickel is controlled by genes with limited effects.

ACK N OWLED G M ENTS
We would like to thank the Natural Sciences and Engineering Research Council of Canada (NSERC) for financial support and Seq Matic (Fremont California) for RNA sequencing. We are also grateful to the National Forest Service, Tree Seed Centre, Canadian Forest Services (Fredericton) for providing certified A. rubrum seeds used in this study.

AVA I L A B I LIT Y O F DATA A N D M ATE R I A L S
We have deposited data in the following depository: