Transcriptome Analysis of Acid-Responsive Genes and Pathways Involved in Polyamine Regulation in Iron Walnut

We reported changes in the co-regulated mRNA expression in iron walnut (Juglans sigillata) in response to soil pH treatments and identified mRNAs specific to acidic soil conditions. Phenotypic and physiological analyses revealed that iron walnut growth was greater for the pH 4–5 and pH 5–6 treatments than for the pH 3–4 and pH 6–7 treatments. A total of 2768 differentially expressed genes were detected and categorized into 12 clusters by Short Time-series Expression Miner (STEM). The 994 low-expression genes in cluster III and 255 high-expression genes in cluster X were classified as acid-responsive genes on the basis of the relationships between phenotype, physiology, and STEM clustering, and the two gene clusters were analyzed by a maximum likelihood (ML) evolutionary tree with the greatest log likelihood values. No prominent sub-clusters occurred in cluster III, but three occurred in cluster X. The Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis indicated that acid-responsive genes were related primarily to arginine biosynthesis and the arginine/proline metabolism pathway, implying that polyamine accumulation may enhance iron walnut acid stress tolerance. Overall, our results revealed 1249 potentially acid-responsive genes in iron walnut, indicating that its response to acid stress involves different pathways and activated genes.


Introduction
Juglans sigillata Dode or iron walnut (whose nut shells have many pits and depressions in their surface) is an important species of walnut tree known for its production of edible nuts and hard wood. This species is sometimes grown in gardens and parks as an ornamental plant. It is commonly found in Yunnan but also in Guizhou, Sichuan, and Tibet in China and grows in soils that have a pH of approximately 4-6 [1]. Due to acid rain pollution and aluminium toxicity [2], the soils in the Long River region and in southern China have gradually become acidified. Acid soil stress affects tree growth and development and endangers the agronomic yields of trees, which have also evolved various molecular and physiological strategies to counteract acid soil stress. Therefore, an in-depth understanding of the gene expression regulatory mechanisms activated in response to acid stress is needed for efficient improvement.
Studies involving acid treatments have been conducted on various species across multiple kingdoms, including Escherichia coli [3], Glycine max [4], Lactococcus lactis [5], Lactobacillus plantarum [6], Listeria monocytogenes [7], and Spinacia oleracea [8]. The ability of plants to mitigate acid rain stress and sustain productivity may be related to the scavenging of stress-induced toxic oxygen species, such as hydrogen peroxide, hydroxyl radicals and superoxide radicals [9]. Catalases and peroxidases such as hydrogen peroxide, hydroxyl radicals and superoxide radicals [9]. Catalases and peroxidases are two major enzyme groups whose members participate in the removal of hydrogen peroxide in plants [10,11]. Moreover, polyamines are well known for their anti-senescence and anti-stress effects because of their antioxidant properties as well as their acid-neutralizing, membrane-stabilizing, and cell-wall-stabilizing abilities [12][13][14][15][16][17]. Polyamines are maintained at optimal levels by antizymes and antizyme inhibitors and perform various cellular functions [18], but at high levels, these molecules are cytotoxic [19]. The stress response (including acid stress) as a potential mechanism may inhibit the proliferation of polyamine-depleted cells [20]. However, the molecular mechanism of how polyamines act in these processes is unclear, and plant acid tolerance remains elusive.
Little research is currently available on tree responses to acid stress. Studies involving iron walnut have focused mainly on the inhibitory effects of the extracts of its green husks [21], the chemical compositions and antioxidant capabilities of oils [22], speciation and phylogenetic relationships [23], its rapid Pleistocene diversification [24], the phenylpropanoid pathway [25], and anthocyanin biosynthesis [26]. Breeding for acid tolerance is particularly challenging because of the genetic complexity of this trait. Iron walnut has a complex genome (2n = 32) [27], ~606 Mb in size [28], and it is well documented that acid tolerance results from cooperative interactions among multiple morphological, physiological, and biochemical characteristics. In this study, we used transcriptome sequencing (transcriptome-seq) to identify differentially expressed genes (DEGs) in response to acid stress; accordingly, we used the leaf tissue of acid-resistant iron walnut landrace plants grown in different soil pH treatments. This analysis serves as a reference for future studies on iron walnut response to various stresses, such as drought, cold, and salt stresses.

Plant Growth under Acidic Soil Conditions
An endemic iron walnut variety that is relatively acid tolerant was used for the analysis of acid tolerance-related genes. The pH of the original soil sample used for planting was measured and adjusted to generate a series of acidic soils (pH 3.25, pH 4.59, pH 5.26, pH 6.54) by adding the appropriate amount of calcium oxide (CaO). Twelve iron walnut seedlings were grown under the same conditions in a growth chamber (14 h photoperiod, 22 °C temperature, and 70% relative humidity; the soil pH was maintained at pH 3-4, pH 4-5, pH 5-6, or pH 6-7 during the acid stress treatment). One seedling was planted per pot, and three pots per soil pH treatment constituted one biological replicate. Seedlings approximately one year old were subjected to acid soil stress treatment for two months ( Figure 1). Afterward, the leaf samples were first used to measure the net photosynthetic rate (Pn) on the basis of photosynthetically active radiation (PAR). The leaves were then harvested, flash frozen in liquid nitrogen and immediately stored at −80 °C for photosynthetic pigment and RNA sequencing (RNA-seq) analyses.

Contents of Photosynthetic Pigments and Pn-PAR Response Curves
To determine the effect of each pH treatment on photosynthesis, Pn-PAR response curves were generated in accordance with the non-rectangular hyperbolic photosynthetic model described by Farquhar et al. [29]. The chlorophyll and carotenoid contents were measured using the acetone-alcohol mixed extraction method described by Romero and Sanchez [30].

RNA Isolation and Iron Walnut RNA-Seq
RNA was extracted with an RNAprep Pure Plant Kit (Polysaccharides & Polyphenolics-rich) (Tiangen Biotech Co., Ltd., Beijing, China) according to the instructions provided. After extracting the total RNA and applying a DNase I treatment, we used magnetic beads with Oligo (dT) to isolate the mRNA. The mRNA molecules were fragmented at 150 bp, and cDNAs were synthesized via mRNA fragments as templates. The cDNA fragments were purified and resolved with EB buffer for end repair, single-nucleotide adenine (A) addition and adapter connections. After PCR amplification, an Agilent 2100 Bioanalyzer and an ABI Step One Plus Real-Time PCR System were used to quantify and qualify the sample library, respectively. The 150 bp library was then sequenced via an Illumina HiSeq™ 4000. Clean reads were obtained using Soap nuke software (http://www.seq500.com/uploadfile/SOAPnuke.zip) by removing "dirty" reads containing adapter sequences, sequences with more than 10% unknown bases ("N"), and low-quality reads (reads with >50% of bases identified as low quality, where low-quality bases were defined as those whose sequencing quality was no more than 5). All of the clean reads were mapped against the genome data of walnut (ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/001/411/ 555/GCF_001411555.1_wgs.5d) using HISAT2 software (version 2.1.0).

DEG Identification and Functional Annotation
Gene expression was calculated by the fragments per kilobase of exon model per million mapped reads (FPKM) method [31]. HTSeq software (version 0.10.0) was used to calculate the gene expression levels via a union model. The genes with FPKM values >1 were subjected to subsequent analysis. The DESeq method [32] was then used to identify DEGs based on a negative quadratic term distribution model. The raw counts and FPKM values were output by HTSeq, and the raw count of each gene was input into DESeq for differential expression analysis. The q-value of the DEGs was less than 0.05 (|Log 2 (FC)| ≥ 1, FC = fold change). R software (version 3.5.3) was used to construct and analyze Venn diagrams, as well as to analyze expression patterns and to perform a cluster analysis with the heatmap clustering method. The DEGs were ultimately annotated by the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases via Omic Share tools (http://www.omicshare.com/).

Acid-Responsive Gene Identification and Enrichment Analysis
We compared the DEG profiles of plants under acid stress via STEM clustering [33]. Significant DEG profiles (p < 0.05) were further analyzed for acid-responsive genes based on the relationship between these gene expression profiles and those of iron walnut seedlings subjected to different soil pH treatments. Molecular phylogenetic analysis of these acid-response-related genes was also performed with MEGA7 [34]; the sequences were aligned by ClustalW, and molecular phylogenetic relationships were inferred by the maximum likelihood (ML) method based on the Tamura-Nei model [35]. Initial trees for the heuristic search were obtained automatically by applying the neighbor-joining and Bio NJ algorithms to a matrix of pairwise distances that were estimated using the maximum composite likelihood approach and then selecting the topology with the superior log likelihood value. The tree is drawn to scale, and the branch lengths represent the number of substitutions per site. The acid-responsive genes related to different biological functions were classified and grouped via GO analysis. KEGG pathway analysis was also used to identify biological pathways of the acid-responsive genes.

Quantitative Real-Time PCR (qRT-PCR) Validation
To validate the results of the RNA-seq data, qRT-PCR was performed to analyze the expression of genes with SYBR Green I using the iron walnut 18S gene as a standard control [36]. cDNAs were reverse-transcribed using a True Script First Stand cDNA Synthesis Kit (Aidlab, Beijing, China). Reverse transcription was performed in a total volume of 20 µL, which consisted of 1 µL of RNA, 4 µL of 5× RT Reaction Mix, 0.8 µL of Rondam primer/oligo (dT), 0.8 µL of True Script H -RTase/RI Mix, and 13.4 µL of RNase-free dH 2 O. The reaction tubes were maintained at 42 • C for 40 min and then at 65 • C for 10 min. The cDNAs were then used as templates for real-time PCR using gene-specific primers. Real-time PCR was conducted in a total volume of 10 µL, which consisted of 5 µL of 2× SYBR ® Green Supermix, 1 µL of primers, 1 µL of cDNA, and 3 µL of ddH 2 O. The reaction conditions were 95 • C for 3 min, followed by 39 cycles of 95 • C for 10 s and 57-61 • C for 30 s. The relative quantification of gene expression was determined by qRT-PCR software 3.0 using the Pfaffl method. Mean values and standard errors were calculated from three independent experiments consisting of three biological replicates of iron walnut, and the data were normalized to the relative efficiency of each primer pair.

Determination of Pn-PAR Response Curves and Photosynthetic Pigment Contents
To investigate the effects of each pH iron walnut in such pH treatment, Pn-PAR response curves were generated (Figure 2a). The trend of the Pn-PAR response curves for the iron walnut seedlings was as follows: pH 4-5 > pH 5-6 > pH 6-7 > pH 3-4. The photosynthetic pigment contents were also measured ( Figure 2b). The total contents of chlorophyll a and b were 14.29%, 11.08%, and 9.09% greater for the iron walnut seedlings at pH 4-5, pH 5-6, and pH 6-7 than at pH 3-4, respectively. Similarly, the content of chlorophyll b was 55.45%, 46.45%, and 25.98% greater for the iron walnut seedlings at pH 4-5, pH 5-6, and pH 6-7, respectively, than at pH 3-4. Moreover, the content of chlorophyll a was 6.02%, 3.97%, and 5.70% greater for the iron walnut seedlings at pH 4-5, pH 5-6, and pH 6-7 than at pH 3-4, respectively. However, the carotenoid content showed the opposite trend from the total chlorophyll content; it was 13.25% and 13.34% lower for the iron walnut seedlings at pH 4-5 and pH 5-6, respectively, but 6.89% greater in the iron walnut seedlings at pH 6-7. In total, the physiological data showed that the growth of iron walnut was greater in the pH 4-5 and pH 5-6 treatments than in the pH 3-4 and pH 6-7 treatments. To validate the results of the RNA-seq data, qRT-PCR was performed to analyze the expression of genes with SYBR Green I using the iron walnut 18S gene as a standard control [36]. cDNAs were reverse-transcribed using a True Script First Stand cDNA Synthesis Kit (Aidlab, Beijing, China). Reverse transcription was performed in a total volume of 20 μL, which consisted of 1 μL of RNA, 4 μL of 5× RT Reaction Mix, 0.8 μL of Rondam primer/oligo (dT), 0.8 μL of True Script H -RTase/RI Mix, and 13.4 μL of RNase-free dH2O. The reaction tubes were maintained at 42 °C for 40 min and then at 65 °C for 10 min. The cDNAs were then used as templates for real-time PCR using gene-specific primers. Real-time PCR was conducted in a total volume of 10 μL, which consisted of 5 μL of 2× SYBR ® Green Supermix, 1 μL of primers, 1 μL of cDNA, and 3 μL of ddH2O. The reaction conditions were 95 °C for 3 min, followed by 39 cycles of 95 °C for 10 s and 57-61 °C for 30 s. The relative quantification of gene expression was determined by qRT-PCR software 3.0 using the Pfaffl method. Mean values and standard errors were calculated from three independent experiments consisting of three biological replicates of iron walnut, and the data were normalized to the relative efficiency of each primer pair.

Determination of Pn-PAR Response Curves and Photosynthetic Pigment Contents
To investigate the effects of each pH iron walnut in such pH treatment, Pn-PAR response curves were generated (Figure 2a). The trend of the Pn-PAR response curves for the iron walnut seedlings was as follows: pH 4-5 > pH 5-6 > pH 6-7 > pH 3-4. The photosynthetic pigment contents were also measured ( Figure 2b). The total contents of chlorophyll a and b were 14.29%, 11.08%, and 9.09% greater for the iron walnut seedlings at pH 4-5, pH 5-6, and pH 6-7 than at pH 3-4, respectively. Similarly, the content of chlorophyll b was 55.45%, 46.45%, and 25.98% greater for the iron walnut seedlings at pH 4-5, pH 5-6, and pH 6-7, respectively, than at pH 3-4. Moreover, the content of chlorophyll a was 6.02%, 3.97%, and 5.70% greater for the iron walnut seedlings at pH 4-5, pH 5-6, and pH 6-7 than at pH 3-4, respectively. However, the carotenoid content showed the opposite trend from the total chlorophyll content; it was 13.25% and 13.34% lower for the iron walnut seedlings at pH 4-5 and pH 5-6, respectively, but 6.89% greater in the iron walnut seedlings at pH 6-7. In total, the physiological data showed that the growth of iron walnut was greater in the pH 4-5 and pH 5-6 treatments than in the pH 3-4 and pH 6-7 treatments.

Illumina Sequencing of Different cDNA Libraries
cDNA libraries were constructed from leaves harvested two months after acid stress; three biological replicates were sequenced using the Illumina HiSeq™ 4000 platform. After cleaning and checking the read quality, we obtained nearly 32.87 million 150 bp paired-end clean reads. Among them, all of the Q30 (Phred score) values were greater than 94%, and the GC content was greater than 45% in each sample. The statistical data of the sequenced samples are shown in Table 1. The clean reads of each sample were mapped to the full gene set of walnut. The transcriptome data from 12 samples mapped to 27,316 genes, which covered 74.11% of the whole gene set. The major mapping reads indicated reliable transcriptome data.

DEG Analysis and Validation of the Iron Walnut Sequencing Data
The gene expression in each sample is presented in Table S1. After removing redundant DEGs between six comparison groups, a total of 2768 DEGs (q < 0.05) were identified based on a pairwise comparison analysis (Figure 3a). Six cross-comparisons between the different gene sets were illustrated in a Venn diagram ( Figure 3b). A cross-comparison of pH 3-4 vs. pH 4-5 showed that the most-expressed genes (748) were unique to the other cross-comparisons, while pH 3-4 vs. pH 6-7 resulted in fewer uniquely expressed genes (2). The heatmap clustering showed the gene expression patterns for the 2768 DEGs to be clustered into three categories ( Figure 3c). The DEGs in the pH 3-4 treatment differed from those in the other three treatments (pH 4-5, pH 5-6, pH 6-7), and the DEGs in the pH 6-7 treatment also differed from those in the other two treatments (pH 4-5, pH 5-6). Among these DEGs, ten were confirmed via qRT-PCR to be distinctly differentially expressed between treatments (Table 2, Figure 4). The results of this experiment were essentially consistent with the RNA-seq data.

GO-and KEGG-Annotated DEGs
The top ten GO term classifications of the DEGs from six cross-comparisons between the different gene sets are illustrated in Figure 5. Among these GO terms, 45% of DEGs belong to the biological process category, 47% of DEGs belong to the molecular function category, and 8% of DEGs belong to the cellular component category. In fact, 40% of the DEGs were associated with electron transport (81 DEGs) and oxidoreductase activity (80 DEGs), 6% of DEGs were associated with 3deoxy-7-phosphoheptulonate synthase activity (3 DEGs), plastid stroma (3 DEGs), chloroplast stroma (3 DEGs), the chloroplast ribulose bisphosphate carboxylase complex (3 DEGs), the ribulose bisphosphate carboxylase complex (3 DEGs), the ribonucleoside-diphosphate reductase complex (2 DEGs), chloroplast parts (3 DEGs), and plastid parts (3 DEGs). Thirty-eight KEGG pathways (false discovery rate (FDR) < 0.05) from six cross-comparisons between the different gene sets are listed in Table 3. Approximately 85% of the DEGs were assigned to metabolism, and the remaining 15% of the DEGs were nearly equally assigned to genetic information processing, environmental information processing, and organismal systems. In fact, 25% of the DEGs focused on metabolic pathways (229 DEGs), and 1% of the DEGs were associated with the proteasome (1 DEG), other glycan degradation (1 DEG), the mRNA surveillance pathway (1 DEG), the biosynthesis of unsaturated fatty acids (1 DEG), and isoflavonoid biosynthesis (1 DEG). Thirty-eight KEGG pathways (false discovery rate (FDR) < 0.05) from six cross-comparisons between the different gene sets are listed in Table 3. Approximately 85% of the DEGs were assigned to metabolism, and the remaining 15% of the DEGs were nearly equally assigned to genetic information processing, environmental information processing, and organismal systems. In fact, 25% of the DEGs focused on metabolic pathways (229 DEGs), and 1% of the DEGs were associated with the proteasome (1 DEG), other glycan degradation (1 DEG), the mRNA surveillance pathway (1 DEG), the biosynthesis of unsaturated fatty acids (1 DEG), and isoflavonoid biosynthesis (1 DEG).

STEM Cluster Analysis of DEGs and Identification of Acid Resistance-Related Genes
STEM cluster analysis was used to estimate all of the DEGs, as shown in Figure 6 and Table S2. A total of 2768 DEGs were filtered, and 2577 DEGs ultimately formed 12 clusters. Four clusters were significantly related to pH treatment (clusters III, VI, X, XII). Furthermore, all 994 genes in cluster III were clearly expressed at low levels in the pH 4-5 and pH 5-6 treatments, while all 255 genes in cluster X were clearly expressed at high levels in the pH 4-5 and pH 5-6 treatments. In these experiments, iron walnut plants grew better in the pH 4-5 and pH 5-6 treatments than in the pH 3-4 and pH 6-7 treatments. Moreover, the physiological data showed that the growth of iron walnut was greater in the pH 4-5 and pH 5-6 treatments than in the pH 3-4 and pH 6-7 treatments. Hence, the 1249 genes in cluster III and cluster X were found to be potential acid-responsive genes.

Molecular Phylogenetic Analysis of Acid Resistance-Related Genes
Evolutionary analyses of the 994 DEGs in cluster III and the 255 DEGs in cluster X were conducted with MEGA7 ( Figure 7). The trees with the highest log likelihood values are shown in Figure 7a,b (−2,370,992.00 and −120,896.32). The 994 low-expression DEGs in cluster III included 1339 nucleotide sequences, as shown in Figure 7a, and the 255 high-expression DEGs in cluster X included 423 nucleotide sequences, as shown in Figure 7b. The total length of the combined sequences in cluster III is 7840 bp and that in cluster X is 5416 bp after alignment by ClustalW in Mega 7. All positions with gaps and missing data were eliminated, then we obtained 1017 bp in cluster III, and 1964 bp in cluster X. There was a total of 208 positions (Figure 7a) and 258 positions (Figure 7b) in the final dataset. No prominent sub-clusters were observed in cluster III (Figure 7a), but three notable sub-clusters were observed in cluster X (Figure 7b). Two long nucleotide sequences with nearly the same number of nucleotides were observed in the green (193 sequences) and red (194 sequences) subclusters. Thirty-six nucleotide sequences formed the small blue sub-cluster. Interestingly, each subcluster in cluster X contained sequences related to photosynthesis and chlorophyll (Figure 7b). In total, differences occurred within high-expression DEGs in cluster X but not within low-expression DEGs in cluster III.

Molecular Phylogenetic Analysis of Acid Resistance-Related Genes
Evolutionary analyses of the 994 DEGs in cluster III and the 255 DEGs in cluster X were conducted with MEGA7 (Figure 7). The trees with the highest log likelihood values are shown in Figure 7a No prominent sub-clusters were observed in cluster III (Figure 7a), but three notable sub-clusters were observed in cluster X (Figure 7b). Two long nucleotide sequences with nearly the same number of nucleotides were observed in the green (193 sequences) and red (194 sequences) sub-clusters. Thirty-six nucleotide sequences formed the small blue sub-cluster. Interestingly, each sub-cluster in cluster X contained sequences related to photosynthesis and chlorophyll (Figure 7b). In total, differences occurred within high-expression DEGs in cluster X but not within low-expression DEGs in cluster III. Figure 7. Evolutionary analyses of cluster III (a) and cluster X (b) conducted with MEGA7. The evolutionary history was inferred using the maximum likelihood (ML) method based on the Tamura-Nei model. Initial tree(s) for the heuristic search were obtained automatically by applying the neighbor-joining and Bio NJ algorithms to a matrix of pairwise distances that were estimated using the maximum composite likelihood approach and then selecting the topology with the superior log likelihood value. The tree is drawn to scale, and the branch lengths represent the number of substitutions per site. All positions containing gaps and missing data were eliminated. The red color represents the cluster including 194 sequences, while the green color represents another cluster including 193 sequences, the blue color represents the rest cluster including 36 sequences.

Gene Function and Pathway Enrichment Analyses of Acid Resistance-Related Genes
The combination of the above strategies resulted in final totals of 994 and 255 DEGs. We then performed GO and KEGG enrichment analyses of these two groups of genes ( Figure 8). The assigned GO functions of the genes in cluster III and cluster X broadly encompassed 45 and 37 functional categories, respectively. The 994 low-expression DEGs in cluster III prominently included the GO terms cell periphery, plasma membrane, oxidoreductase activity, transporter activity, anion transport, oxidation-reduction process, transmembrane transport, and defence response (Figure 8a). The 255 high-expression DEGs in cluster X were represented by the GO terms thylakoid, chloroplast, plastid, photosystem I, envelope, extracellular region, photosynthesis, light reaction, organic acid biosynthetic process, and carboxylic acid biosynthetic process (Figure 8b). Approximately half of the GO functional enrichment terms overlapped between clusters III and X, including response to stimulus, binding, and extracellular region. Furthermore, oligopeptide transport, peptide transport, chitin binding, oxidoreductase activity, chloroplast stroma, and cell cortex were involved in the scavenging of stress-induced toxic oxygen species and polyamine regulation.
The KEGG pathway analysis revealed that the 994 low-expression DEGs in cluster III were involved mainly in eight pathways: Limonene/pinene degradation, stilbenoid/diarylheptanoid/ gingerol biosynthesis, galactose metabolism, selenocompound metabolism, anthocyanin biosynthesis, glyoxylate/dicarboxylate metabolism, protein processing in the endoplasmic reticulum, and plant-pathogen interactions (Figure 8c). The analysis also revealed that the 255 high-expression DEGs in cluster X were involved mainly in another eight pathways: Photosynthesis, fatty acid degradation, α-linolenic acid metabolism, biosynthesis of amino acids, alanine/aspartate/glutamate metabolism, glycine/serine/threonine metabolism, arginine biosynthesis, and arginine/proline metabolism (Figure 8d). There were no overlapping KEGG pathways between clusters III and X. Among these pathways, glyoxylate/dicarboxylate metabolism, biosynthesis of amino acids, alanine/aspartate/glutamate metabolism, glycine/serine/threonine metabolism, arginine biosynthesis, and arginine/proline metabolism were involved in polyamine metabolism and regulation.

Discussion
Acidic soil is prevalent in the Yangtze valley and southern area of China. Among the major trees of economic value, iron walnut is relatively tolerant to acid stress and therefore is often grown on more marginal sites. As a landrace, iron walnut has developed a strong tolerance and adaptability to acid stress. Using high-throughput RNA-seq technology, we compared in detail the transcriptional differences and overlap in response to different levels of soil pH and identified the changes in gene expression in the iron walnut seedlings during the whole acid response process.

Acid Resistance-Related Gene Identification
In this study, the mRNA of iron walnut seedlings exhibiting good acid tolerance was sequenced using an Illumina HiSeq™ 4000 platform with Sera-Mag Magnetic Oligo-dT Beads. A clear bioinformatic map of the mRNAs involved in multiple biological processes was produced. As a result, 32.87 million clean reads were collected from 12 samples of plants subjected to different soil pH treatments; this number met the requirements for further analysis. The clean Q30 reads constituted more than 94% of the total, suggesting high-quality sequencing. HISAT2 software was used to BLAST the transcriptome data against the walnut genome; 74% of the reads were mapped to the reference genome. The major mapping reads indicated reliable transcriptome data. These non-mapped tags most likely represent regions where the reference genome is incomplete [37]. Alternatively, there are genetic differences between the reference genome and the iron walnut seedlings used in this study [28]. Moreover, the RNA-seq data for reference genome annotations might not represent all major tissue types, developmental stages, or statuses during responses to abiotic and biotic stresses [38].
The global analysis of the DEGs provided a comprehensive dataset on the response to acid soil stress in the leaves of iron walnut seedlings. We identified 12 clusters for all DEGs and obtained four clusters that were significantly related to soil pH (p < 0.05). The genes in cluster III and cluster X were determined to be potential acid-responsive genes on the basis of the relationship between these gene expression profiles (STEM clustering trends) and the growth of the iron walnut seedlings subjected to the pH 4-5 and pH 5-6 treatments (the phenotype), as well as the Pn-PAR response curves and photosynthetic pigment contents (physiological data). Similarly, potential genes in both maize and humans have also been identified by gene expression profiling and phenotyping [39][40][41].

Analysis of Acid Resistance-Related Genes
The acid resistance-related genes in clusters III and X were further analyzed by the ML method. No prominent sub-clusters were observed in cluster III, but three notable sub-clusters were observed in cluster X, indicating a large difference between the low-expression genes in cluster III and the high-expression genes in cluster X. Two nearly identical large sub-clusters and a small sub-cluster in cluster X all contained sequences related to photosynthesis and chlorophyll, indicating that differences in nucleotide sequences occurred within the high-expression genes in cluster X, although these nucleotide sequences have similar biological functions. A similar phenomenon was also found in Arabidopsis thaliana [42], Capsicum annuum [43], Lactobacillus kefiranofaciens [44], Rana arvalis [45], and Zea mays [46].
The GO functional terms, including plasma membrane, oxidoreductase activity, anion transport, oxidation-reduction process, transmembrane transport, and defence response, were the most affected among the terms associated with the low-expression genes in cluster III, indicating that the genes in these different functional categories were predicted to be involved at low expression levels in iron walnut seedlings. GO terms including chloroplast, photosystem I, photosynthesis, and light reaction were associated with the most strongly affected, high enough genes in cluster X, indicating that the genes associated with these different functional categories can be predicted to be involved at high expression levels in iron walnut seedlings. GO terms such as photosystem I and photosynthesis were reported by Shikanai et al. [47], Guan et al. [48], and Otani et al. [49]. Approximately half of the GO functional enrichment terms overlapped between clusters III and X, including response to stimulus, binding, and extracellular region, suggesting that these associated genes can be predicted to be involved in essential plant biological processes [17,50,51]. Several similar GO terms, such as oligopeptide transport, peptide transport, chitin binding, oxidoreductase activity, chloroplast stroma, and cell cortex, have been identified as being associated with anti-senescence and anti-acid stress effects in A. thaliana, C. annuum, Nicotiana tabacum, Triticum aestivum, and Z. mays [11,12,16,40,46,52], implying that these GO terms are involved in the scavenging of stress-induced toxic oxygen species and polyamine regulation.
KEGG pathway analysis in this study revealed sixteen pathways that were significantly affected by acid soil stress. The low-expression genes associated with eight pathways, glyoxylate/ dicarboxylate metabolism, limonene/pinene degradation, plant-pathogen interaction, etc., in cluster III significantly differed throughout the whole acid soil stress period, implying that these genes are involved in functions of defence, antioxidation, free radical scavenging, sugar, lipids as an energy source, protein processing, and polyamine metabolism and regulation. The expression of high-expression genes associated with another eight pathways (alanine/aspartate/glutamate metabolism, arginine biosynthesis, arginine/ proline metabolism, etc.) in cluster X significantly varied throughout the whole acid soil stress period.
These results indicated that polyamine accumulation enhanced iron walnut acid stress tolerance. Similar studies by Gupta et al. [11] reported that the accumulation of polyamines occurs under many types of abiotic stress, including acid stress, and modulation of their biosynthetic pathway results in tolerance to stress. Hussain et al. [9] also reported a strong correlation between the presence of polyamines and improved tolerance to environmental stresses (acid drought, salinity, low temperature, oxidative stress and metal toxicity). Taie et al. [16] reported that polyamines promote growth and antioxidant activity and protect genomic DNA in heavy metal-stressed wheat plants. Similarly, Liu et al. [53] noted that polyamines provided environmental stress tolerance to plants. Majumdar et al. [46] demonstrated that polyamines contributed to the resistance of maize against Aspergillus flavus infection and to aflatoxin production. Amrani et al. [15] revealed that polyamines improved sucrose-induced tolerance to atrazine-mediated chemical stress in A. thaliana. Furthermore, Velikova et al. [4] proposed that exogenous polyamines could be attributed to protecting acid rain-treated bean plants via acid-neutralizing and antioxidant effects, as well as by their ability to stabilize membranes by associating with negatively charged phospholipids.
The hypothetical regulatory mechanisms of enhanced acid stress tolerance via polyamines are illustrated in Figure 9. These varied hypothetical mechanisms were proposed by Liu et al. [53], Hussain et al. [9], Gupta et al. [11], Taie et al. [16], and Zarza et al. [17]. The biosynthesis of polyamines in plants has been well demonstrated [59][60][61][62][63][64][65]. Putrescine is produced either directly from ornithine or indirectly from arginine. Putrescine is converted into spermidine, and spermidine is converted into spermine, both of which require methionine. Polyamines include putrescine (diamine), spermidine (triamine), and spermine (tetramine). The hypothetical roles of polyamines include antioxidant responses, programmed cell death, acting as a compatible solute in higher plants, acting as a signaling molecule during stress responses, regulating ion channels and Ca 2+ homeostasis, and other unknown pathways. Previous works [4,9,11,[15][16][17][18]20,43,46,52,53,62] have indicated that accumulated polyamines may play vital roles in counteracting stress (1) as antioxidants, free radical scavengers, and membrane stabilizers that can interact freely with DNA/RNA/proteins/membrane lipids to prevent the membrane system from denaturing under acid stress conditions; (2) through the cationic nature of spermine +4 > spermidine +3 > putrescine +2 at physiological pH and/or improvement of the K + /Ca 2+ ion balance; (3) as cellular signals in intrinsic cross-talk with abscisic acid (ABA), hydrogen peroxide, and nitrous oxide during acid stress responses by a complex network; and (4) by increasing cell proliferation, decreasing apoptosis and increasing the expression of genes associated with tumor invasion and metastasis, thus making their metabolism a target for cancer treatment and prevention.
Genes 2019, 10, x FOR PEER REVIEW 17 of 21 Figure 9. Overview of hypothetical regulatory mechanisms of enhanced acid stress tolerance via polyamines. Acid stress induces the accumulation of polyamines, which then function in antioxidant responses, in programmed cell death, as compatible solutes in higher plants, as signaling molecules during stress responses, in the regulation of ion channels and Ca 2+ homeostasis, and in other unknown pathways, by scavenging excess free radicals or stabilizing membranes by binding to negatively charged groups, leading to enhanced acid stress tolerance.

Conclusions
In summary, this study provides a comprehensive analysis of the acid-responsive genes and transcriptome expression profiles of iron walnut leaves via a combination of RNA-seq, STEM profiling, the ML method, and physiological approaches. Our results revealed 1249 potentially acidresponsive genes in iron walnut, which we annotated by using the GO and KEGG databases, showing that polyamine accumulation enhanced iron walnut acid stress tolerance. These results fill the knowledge gap in the available literature concerning acid resistance-related genes in iron walnut. Next, we seek to understand how environmental changes (such as acid stress) are communicated via cellular signal transduction to induce a coordinated metabolic response and how the function of polyamines is influenced by transcriptional and post-translational modifications.

Supplementary Materials:
The following are available online at www.mdpi.com/xxx/s1, Table S1: Gene expression for each sample. The first column in the table presents the   Acknowledgments: The authors thank Qianwen Xiao and Guanglan Pu for their assistance in obtaining plant materials and support. Mauritius

Conflicts of Interest:
The authors declare that they have no conflicts of interest. Figure 9. Overview of hypothetical regulatory mechanisms of enhanced acid stress tolerance via polyamines. Acid stress induces the accumulation of polyamines, which then function in antioxidant responses, in programmed cell death, as compatible solutes in higher plants, as signaling molecules during stress responses, in the regulation of ion channels and Ca 2+ homeostasis, and in other unknown pathways, by scavenging excess free radicals or stabilizing membranes by binding to negatively charged groups, leading to enhanced acid stress tolerance.

Conclusions
In summary, this study provides a comprehensive analysis of the acid-responsive genes and transcriptome expression profiles of iron walnut leaves via a combination of RNA-seq, STEM profiling, the ML method, and physiological approaches. Our results revealed 1249 potentially acid-responsive genes in iron walnut, which we annotated by using the GO and KEGG databases, showing that polyamine accumulation enhanced iron walnut acid stress tolerance. These results fill the knowledge gap in the available literature concerning acid resistance-related genes in iron walnut. Next, we seek to understand how environmental changes (such as acid stress) are communicated via cellular signal transduction to induce a coordinated metabolic response and how the function of polyamines is influenced by transcriptional and post-translational modifications.  Acknowledgments: The authors thank Qianwen Xiao and Guanglan Pu for their assistance in obtaining plant materials and support.

Conflicts of Interest:
The authors declare that they have no conflicts of interest.