The maize (Zea mays L. inbred line Zheng 58) seeds were surface sterilized, washed with ddH₂O, and then germinated overnight in an incubator at 30°C. The seedlings were planted in a growth chamber with a photoperiod of 16-h light/8-h dark and a relative humidity of 60%. A half-strength, modified Hoagland nutrient solution (pH = 5.8) was used and changed every 3 days. The seedlings were grown in two groups of 40 × 40 cm pots (10 seedlings per pot) representing different treatments. Two-week-old seedlings were grown in nutrient solution with or without 100 μM cadmium chloride for 7 days. Then the root samples were harvested and stored at -80°C in preparation for further assays. Total RNA was extracted from the different samples using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s protocol. RNA contamination was detected using a Bioanalyzer 2100 and RNA 6000 Nano LabChip Kit (Agilent, Santa Clara, CA, USA) and removed by 1% agarose gel electrophoresis.

Approximately 10 μg of total RNA was subjected to isolate poly (A) mRNA with poly-T oligo attached magnetic beads (Invitrogen, Beijing, China). Following purification, the mRNA was fragmented into small pieces using divalent cations at elevated temperatures. Then the cleaved RNA fragments were reverse-transcribed to create the final cDNA library in accordance with the protocol for the mRNA-Seq sample preparation kit (Illumina, San Diego, CA, USA). The average insert size for the paired-end libraries was 300 bp (±50 bp). Then we performed paired-end sequencing on an Illumina Hiseq2000/2500 (LC Sciences, USA) following the manufacturer’s protocol. Six RNA libraries consisted of three control libraries and three Cd-treated libraries.

We used the Illumina paired-end RNA-seq approach to sequence the maize root transcriptome, which generated 244 million paired-end reads. Three repetitions of each treatment have been sequenced. This yielded 30.37 Gb of sequences, which was approximately 13.2 times the size of the genome (2.3 Gb) (Schnable et al., 2009). Prior to assembly, low quality reads, including reads containing sequencing adaptors, sequencing primers, and nucleotides with a q quality score lower than 20, were removed. The raw sequence data have been submitted to the NCBI Short Read Archive with an accession number of GSE74516.

We aligned the reads of the different samples to the MaizeGDB¹ maize reference genome using Tophat package v2.0.9 (Trapnell et al., 2010), which initially removes a proportion of the reads based on the quality information accompanying each read, and then maps the reads to the reference genome. Tophat allows multiple alignments per read and a maximum of two mismatches when mapping the reads to the reference genome. Tophat builds a database of potential splice junctions and confirms these by comparing the previously unmapped reads against the database for putative junctions.

The aligned read files were processed by Cufflinks, which uses the normalized RNA-seq fragment counts to measure the relative abundances of the transcripts. Cufflink was used to de novo assemble the transcriptome, and then Cuffmerge was used to integrate all the transcripts from the different samples to generate unique transcripts. The final unigenes that showed differential expressions between the different maize treatments were detected by DEGseq software using three replicates per sample (Anders and Huber, 2010). The unigene expression levels were calculated using reads per kilobase per million reads (RPKM), which eliminated the influences of gene length and sequencing level during the calculation of gene expression. A general Chi-squared test of statistical significance was used, and the false discovery rate (FDR) for the results was controlled (FDR < 0.05) (Ran and Peng, 2016). Significantly altered genes were described using heatmap analysis with unsupervised hierarchical clustering. The raw intensity (RPKM) was log₂ transformed and then used to calculate the Z scores (Cheadle et al., 2003).

To assign putative functions to differentially expressed genes (DEGs) during maize responses to the treatments, various bioinformatics approaches were used for further annotation, classification, and metabolic pathway analysis. First, the DEGs were aligned to the web-based agriGO tools. GO enrichment analysis of DEGs was implemented using singular enrichment analysis (SEA) by comparing a query list of DEGs to a background gene set (FDR < 0.05). The SEACOMPARE tool was used for a comparative analysis that integrated the SEA cross-comparison. Finally, KOBAS software was used to test DEG statistical enrichment in the KEGG pathways (Xie et al., 2011).

Several DEGs with putative functions were selected randomly and validated by qRT-PCR to confirm the results of the RNA-Seq. The primers used in the qRT-PCR experiments were designed by Primer5 software and are listed in Supplementary Table S1. The methods, including RNA extraction from the maize seedling root samples, reverse transcription, and qRT-PCR, were performed according to the manufacturer’s protocols (Clontech, Dalian, China). Briefly, 1 μl of a 1/10 dilution of cDNA in double distilled water was added to 5 μl of 2× UltraSYBR. Then 100 nM of each primer was added to water to make up a final volume 10 μl. The procedure for the PCR were as follows: 95°C for 10 min; 40 cycles of 95°C for 15 s, and 60°C for 60 s. Heat map representation was performed using the average Ct value, and ClustalW software and Treeview were used to visualize the qRT-PCR analysis data. Five biological repeats were used for the expression analyses and the values shown in the figures represent the average values of these five repeats.

Five seedlings for each treatment were oven-dried at 80°C for 3 days and then digested with nitric-perchloric acid (3:1, v/v) at 100°C. A microwave oven (LWY-84B Shenglan, Jiangsu, China) equipped with an infrared ray generation device was used to eliminate any interfering organic substances. The Cd content measurements were made by a graphite furnace atomic absorption spectrophotometer (ICE^TM 3300 AAS, ThermoFisher, USA).

After 7 days of treatment, the dithizone method (Balestri et al., 2014) was used to histochemically detect the Cd in the maize roots. The reddish colored precipitates produced by the dithizone-Cd reaction were determined and analyzed, and the results were used to localize the heavy metal Cd in the maize roots under the different treatments. Five independent roots were collected from each treatment, and the samples were stained for 1.5 h with a dithizone working solution (30 mg dissolved in 60 ml acetone and 20 ml distilled water), washed in clean water, and immediately analyzed using a Carl Zeiss LSM510 laser scanning system.

For hormone treatments, 2-week-old seedlings were transferred to hormone free nutrient solution or nutrient solution with 0.1 μM 1-Naphthaleneacetic acid (NAA) and 1 μM 1-Naphthoxyacetic acid (1-NOA) for 7 days respectively. For Cd treatment, 2-week-old seedlings were grown in nutrient solution with 50 μM cadmium chloride for 7 days. For hormone-Cd combined treatments, 2-week-old seedlings were grown in 0.1 μM NAA + 50 μM cadmium chloride for 7 days or 1 μM 1-NOA + 50 μM cadmium chloride for 7 days. Seedlings in hormone free and Cd free nutrient solution were used as control.

The root samples from control and Cd-treated seedlings were cut and homogenized by 50 mM Tris-HCl buffer, pH 7.6. Then, samples were collected by centrifugation at 12,000 g in a 1.5 ml centrifuge and keep in liquid nitrogen immediately. Five independent biological replicates of 20 mg each were purified after addition of 250 pg of ¹³C6-IAA internal standard using ProElu C18², Auxin content were measured with FOCUS GC-DSQII (Thermo Fisher Scientific Inc., Austin, TX, USA). The IAA oxidase activity was determined according to Xu et al. (2010).

Differences between values were calculated using one-way analysis of ANOVA with Student’s t-test at a significance level of 0.05 in Excel software. All experiments were performed for five biological repeats and the values shown in figures represent the average values of five repeats.

The global gene expression profiles were surveyed using the RNA-Seq approach in order to identify the DEGs that were responsive to Cd stress in the maize seedling roots. The raw Illumina sequencing reads were qualified and adapter trimmed to yield a total of 244 million clean short reads, which contained 30.37 Gb of sequence data from six complementary cDNA libraries. Over 98.76% of the clean reads had quality scores at the Q20 level and over 83% of the clean reads had scores at the Q30 level. A high proportion of the valid, clean reads (53.48–66.24%) were readily mapped onto the maize reference genome sequence in B73 after the different treatments (Supplementary Table S2). The transcriptional abundances of the genes were quantified using Cufflinks and measured as RPKM. A total of 88,557 transcripts have been identified in maize by previous studies, of which, 85,663 transcripts were expressed in our samples.

Transcriptional changes in response to Cd stress were determined by comparing the control and treated maize seedlings transcriptomes using cuffdiff software. The global comparisons of the gene expression profiles for the control and Cd treated samples are shown in Figure 1A, and differences in the expression levels between the two samples are shown as log₂-transformed ratios. After using the 1.5-fold and padj < 0.01 criteria to select genes, 6342 genes, including 3778 induced genes and 2560 down-regulated genes, were identified as being differentially expressed after Cd treatment (Figures 1B,C).

We analyzed the GO terms represented by these genes to obtain useful information about the DEG responses to Cd treatment. In total, 908 enriched GO terms were identified within the DEGs. GO term enrichment analysis indicated that various biological processes and molecular functions, such as DNA binding transcription factor, S-adenosylmethionine biosynthesis, response to oxidative stress, response to biotic stimulus, and peroxidase activity, were significantly enriched in the DEGs (Figure 2; Supplementary Table S3). We mapped the DEGs to the reference canonical pathways in the Kyoto Encyclopedia of Genes and Genomes (KEGG) to further identify the active metabolism pathways involved in the responses to Cd. All DEGs could be classified into 198 predicted biosynthesis pathways, of which, 22 metabolic pathways were significantly enriched (p < 0.05) (Supplementary Table S4). A high proportion of the up-regulated DEGs were enriched in six KEGG pathways. These were plant-pathogen interaction (11%), phenylpropanoid biosynthesis (9%), starch and sucrose metabolism (6%), amino sugar and nucleotide sugar metabolism (6%), cell cycle (5%), and phenylalanine metabolism (4%). In contrast, a large number of down-regulated DEGs were also enriched in the phenylpropanoid biosynthesis (6%), glbenoid, diarylheptanoid and gingerol biosynthesis (6%), glycolysis/gluconegenesis (5%), glutathione metabolism (4%), cysteine and methionine metabolism (4%), and plant-pathogen interaction (4%) KEGG pathways (Figure 3).

The expressions of the auxin signaling pathway genes were analyzed to determine the involvement of auxin and auxin signaling in the maize response to Cd. In our study, a large number of auxin-related genes were identified as DEGs. Comparison of the transcript abundances for auxin biosynthesis, transport, and downstream response genes revealed a universal expression response under Cd treatment (Figure 4A). Five biosynthesis-related genes (ZmYUC2, ZmYUC3, ZmYUC8, ZmYUC9, and ZmYUC10) were significantly down-regulated by Cd treatment; the ZmPIN1 and ZmPIN5 auxin efflux carriers were up-regulated; and ZmPIN4 was down-regulated. The ZmLAX2 and ZmLAX3 auxin influx carriers were induced by Cd treatment; whereas ZmLAX1 was reduced by Cd treatment. Interestingly, only ZmIAA25, ZmARF19, ZmGH3.1, and ZmGH3.9 were up-regulated by Cd treatment, and most of the auxin downstream response genes were down-regulated by Cd treatment.

We performed qRT-PCR assays with independent samples (roots from the control and Cd treated seedlings) to verify the DEGs related to auxin signaling that were identified using RNA-Seq. We randomly selected 12 unigenes from the auxin signaling pathway to validate the RNA-Seq data. The expression levels of these selected genes were basically consistent with the RNA-Seq results (Figure 4B).

We further examined the effects of Cd stress on auxin levels and the endogenous IAA contents in the control and Cd-treated maize seedlings. Compared to the control roots, treatment with 25 μM Cd reduced the endogenous IAA content from 12.6 ng.g^-1 to 5.1 ng.g^-1, and the IAA content was reduced to 4.8 ng.g^-1 by 50 μM Cd (Figure 5A). Furthermore, the relative IAA oxidase activities were measured in the Cd-treated samples and compared to the control samples. The results suggested that IAA oxidase activities were clearly induced by the 25 μM Cd and 50 μM Cd treatments (Figure 5B).

The roles of auxin signaling in root system development have been established by many previous studies (Liu et al., 2009; Wang et al., 2009; Druege et al., 2016). In our study, 0.1 μM NAA and 1 μM 1-NOA (auxin transport inhibitor) was used to examine the roles of auxin and auxin transport in the root system response to Cd stress in maize. The results indicated that an application of 50 μM Cd clearly inhibited the growth of the root system. However, root growth inhibition by 50 μM Cd was markedly alleviated by the addition of 1-NOA, but was increased by the addition of NAA (Figure 6a). The changes in root biomass were also measured under different conditions. These were the control, NAA, 1-NOA, Cd50, NAA + Cd50, and 1-NOA + Cd50 treatments. Low concentrations of NAA and 1-NOA did not significantly change root biomass compared to the control roots. However, 50 μM Cd caused an acute decline in maize root biomass, but this was reversed by the application of 1-NOA and increased by the addition of NAA (Figure 6b).

In situ localization revealed the presence of Cd through the accumulation of reddish precipitates. Few reddish precipitates were observed in the three Cd-untreated roots (control, NAA, and 1-NOA treatments), but large amounts of reddish precipitates accumulated in the three Cd-treated roots (Cd50, Cd50 + NAA, and Cd50 + 1-NOA). More reddish precipitates were observed in the Cd-treated roots than in the control roots. Interestingly, more reddish precipitates were observed in the NAA + Cd50 roots than in the Cd50 roots, which suggested that Cd accumulation was high in the NAA + Cd50 roots. Furthermore, fewer reddish precipitates were observed in the 1-NOA + Cd50 roots than in the Cd50 roots, which indicated that Cd accumulation was low in the 1-NOA + Cd50 roots (Figure 6c). The analytical determination of Cd levels in maize roots under different conditions confirmed the in situ localization of Cd (Figure 6c).