Physiological and Transcriptomic Analyses Reveal Commonalities and Specificities in Wheat in Response to Aluminum and Manganese

1. Introduction

Acid soils are found throughout the world. It is estimated that approximately 30% ~40% of the world’s total land area and more than 50% of the world’s potential arable land have a pH below 5.5; moreover, soil acidification is increasing worldwide. The primary limitations of acid soils are toxic levels of aluminum (Al) and manganese (Mn), as well as suboptimal levels of phosphorous (P).

Aluminum (Al) makes up 7% of the Earth’s crust and is the most abundant metal element. It is fixed in minerals or bound to mineral surfaces in non-acid soils, thus benign to plant growth. When soil acidifies and pH drops below 5, Al³⁺ is solubilized into the soil solution. The activity of Al³⁺ in soil solution increases with the decreasing pH and peaks at pH 4.1 [1]. The primary symptom of Al toxicity is a rapid (beginning within minutes) inhibition of root growth, resulting in a reduced and damaged root system and limited water and mineral nutrient uptake [2].

Manganese is an essential micronutrient and a toxicity element for higher plants, with a narrow concentration window between deficiency and toxicity. With decreasing pH, the amount of exchangeable manganese, mainly in Mn²⁺ form, increases in the soil solution. Mn²⁺ can be readily transported into the root cells and translocated to the shoots, where it is finally accumulated [3]. Mn toxicity symptoms are localized to the shoot, characterized by stunted growth, chlorosis, and necrotic lesions in the leaves, accompanied by the inhibition of photosynthesis and the accumulation of excessive reactive oxygen species (ROS) [2,4].

As the co-existing metal stresses, Al and Mn exert some common toxicity on morphology and physiology, such as inhibiting root growth, disturbing ions balance, and overproducing reactive oxygen species, thus subsequently resulting in oxidative damage, decreasing chlorophyll content, and deteriorating photosynthetic performance [2,5,6]. On the other hand, plants have evolved external exclusion and internal tolerance mechanisms to cope with both active Al and excess Mn. External exclusion mechanism includes the chelation of Al and Mn in the rhizosphere by organic acids and second metabolites [7] and the induction of a pH barrier in the rhizosphere [8,9,10]. By contrast, internal tolerance mechanisms enable plants to deal with Al³⁺ and Mn²⁺ once they enter plant cells, either by forming harmless complexes with organic acids, proteins, phenols, and other organic compounds ligands, by sequestering them to organelles (e.g., vacuoles, endoplasmic reticula, and Golgi bodies), or by rapidly repairing any lesions incurred, including those resulting from oxidative stress [6,9,11].

Due to differences between Al and Mn in the binding site, distribution, and biological functions in plants, some plants’ responses (toxicity and tolerance) to Al and Mn are specific [7,11,12] or similar but play a different role and are dominated by various genes. For example, plants transport Al and Mn from root to leaf or internally across the cytosol into various organelles, apoplast, and cell walls through the activities of series transporters, representing a diverse set of metal transporter families [11,13,14]. In addition, the Complex gene regulatory network of plant hormones, signal transduction, and transcription factors endows plants with specific and common tolerance mechanisms under Al and Mn stress, such as regulating root growth, cell wall modification, and organic acid secretion [11,15,16].

Traditional strategies to maintain crop production on acid soils include the application of lime to raise the soil pH and planting tolerant crop species and cultivars. Although the liming of soil can ameliorate some Al and Mn toxicity and increase calcium and magnesium concentrations in many acid soils worldwide, it is expensive and ineffective in the subsoil. In some cases, heavy lime application may cause a deleterious effect on soil structure and reduce the availability of other cationic micronutrients essential for plant growth [17,18]. Plant species vary in their tolerance to stresses in acid soils. Native plants usually have a more substantial tolerance to adverse factors than crops [19]. Moreover, there is a wide range of tolerance to either Al [20] or Mn [21] among genotypes within the same species. Ideal crop cultivars growing in acid soil should tolerate both Al and Mn toxicity simultaneously. However, Plants may employ some specific mechanisms to tolerate Al or Mn despite their coexistence. Many cultivars tolerate only one toxicity [3,21]. Two strategies could be used to obtain duel-tolerant varieties. One approach is to select and breed manganese-tolerant varieties from aluminum-tolerant germplasms, or vice versa, but it is generally restricted by genetic bottleneck. Another strategy is to design smart plants and pyramid genes to fortify the common and compensate for the differential tolerance. Uncovering the commonalities and the specificities of Al and Mn toxicity and tolerant mechanisms is a prerequisite for identifying and pyramiding hub genes and breeding multi-tolerant genotypes. Although the physiological mechanisms by which plants coordinately adapt to multiple stresses in acid soils have been extensively studied [3,22,23], the molecular mechanisms are still largely unknown. Many genes may have co-evolutionary mechanisms in their regulatory pathways and function in the context of multiple stresses.

In the face of increasing environmental challenges, understanding the molecular and physiological responses to abiotic stresses is paramount for ensuring global food security. Aluminum (Al) and manganese (Mn) stresses pose significant threats to crop productivity [3], particularly in wheat, a staple food for a large portion of the world's population [24]. High-throughput RNA sequencing will be employed to capture the dynamic changes in gene expression profiles, and the transcriptome analysis will provide a genome-wide perspective on the molecular events occurring in wheat plants under Al and Mn stresses [25,26]. This study integrated physiological and transcriptome analysis to reveal the specificity and commonality in the toxicity and tolerance mechanisms of wheat under Al and Mn stresses, providing a new perspective on the molecular network of plants responding to multiple stresses in acid soil and further discovering candidate genes for wheat Al and Mn tolerance, providing genetic resources.

2. Materials and Methods

2.1. Plant Materials and Treatment

Wheat near-isogenic (over 95%) lines differing in Al tolerance at the ALMT1 locus, ET8 (homozygous Al-tolerant) and ES8 (homozygous Al-sensitive) [27], as well as their parents Carazinho (Al tolerant) and Egret (Al sensitive), were used as plant materials. The ET8 and ES8 lines were derived from a cross between Carazinho and Egret, with the resulting progeny backcrossed eight times to Egret or derivatives of Egret. Grain yield increased in Egret-derived lines when Al tolerance from Carazinho was introduced. The seeds were provided by Peter Ryan from CSIRO Plant Industry, Canberra, Australia.

Wheat seeds were surface-sterilized with 0.5% NaClO (v/v) solution for 20 minutes, rinsed with sterile water, and allowed to germinate in wet paper tower for 4~5 days at 4°C before the seedlings were transplanted into continuously-aerated Magnavaca hydroponic solution (pH 4.2) for a 2-day acclimation to low pH in a walk-in growth chamber with a stable temperature of 23°C, a 14 h photoperiod of 400 μmol m⁻² s⁻¹ illumination level, and a relative humidity of 60-80%. Then, the seedlings were stressed by replacing fresh Magnavaca solutions with an addition of 100 µM KAl(SO)₄ (as Al toxicity), 500 µM MnCl₂ (as Mn toxicity), or no addition (as control). Free Al³⁺ and Mn²⁺ concentrations in the solution were 25 µM and 267 µM, respectively, as predicted by the WinIAP program (https://www.sequentix.de/index.php#Winiap). The experiment was arranged in a completely randomized block design with at least three replicates per batch.

2.2. Relative Root Length

Before and four days after stress, the fresh roots of different treatment groups were collected and completely expanded onto the scanner platform (Epson Expression 1000XL). The total length per plant was evaluated by image analysis with WinRhizo software. The RRL, the ratio of net root growth under Al or Mn treatment compared to control, was calculated.

2.3. SPAD and Brown Spots

After four days of stress, the relative content of chlorophyll in fully developed leaves, as dictated by SPAD (Soil and Plant Analyzer Development) value, was determined with a portable Konica SPAD-502 Plus instrument. Brown spots on leaves under Mn stress were checked and photographed.

2.4. Metal Ions Accumulation and Subcellular Distribution

Al and Mn concentrations in roots and shoots and Ca, Mg, Zn, and Fe concentrations in whole-plant were determined after four days of stress. After being dried and pulverized, the plant materials were conventionally digested with mixed acid (concentrated nitric acid: perchloric acid = 4:1) and then determined by inductively coupled plasma atomic emission spectrometry (ICP-AES).

Cell sap in root apices preparation and the ions (Al, Mn) extraction from the residual cell wall were conducted according to [28]. Al and Mn concentrations were determined by ICP-AES.

2.5. RNA Extraction, Library Preparation, Sequencing, and Read Mapping

Total RNA was extracted from 1 cm root apices using TRIzol Reagent (Invitrogen) according to the manufacturer’s instructions. The RNA quality was determined by 5300 Bioanalyser (Agilent) and quantified using the ND-2000 (NanoDrop Technologies). Only high-quality RNA sample (OD260/280 = 1.8~2.2, OD260/230 ≥ 2.0, RIN ≥ 6.5, 28S:18S ≥ 1.0, > 1μg) was used to construct sequencing library and verified by qRT-PCR.

RNA purification, reverse transcription, library construction, and sequencing were performed at Shanghai Majorbio Bio-pharm Biotechnology Co., Ltd. (Shanghai, China) according to the manufacturer’s instructions (Illumina, CA). Thirty-six RNA-seq transcriptome libraries (4 varieties × 3 treatments × 3 biological replicates) were prepared and sequenced.

The raw paired-end reads were trimmed and quality controlled by Fastp [29] with default parameters. The clean reads were separately aligned to the wheat genome with orientation mode using HISAT2 software [30]. The mapped reads of each sample were assembled by StringTie [31] in a reference-based approach. All assembled genes were annotated against Kyoto Encyclopedia of Genes and Genomes (KEGG), Swiss-Prot, Protein family (Pfam), Gene Ontology (GO), Clusters of Orthologous Groups of Proteins (COG), and NCBI non-redundant protein sequences (NR) libraries.

2.6. Differential Expression and Functional Enrichment Analysis.

The expression level of each transcript was calculated according to the fragments per kilobase of exons model per million mapped reads (FPKM) method. RSEM was used to quantify gene abundances [32]. Differential expression analysis was performed using the DESeq2 [33]. DEGs with |log₂FC| ≥ 1 and FDR ≤ 0.001 were considered significantly different expressed genes. Cluster heat diagrams were drawn using Toolkit for Biologists (TB) tools with default settings [34]. Those data were added to the Comprehensive Gene Expression Database with the accession number PRJNA1031207. In addition, GO annotation and KEGG pathway analysis were carried out by Goatools and KOBAS [35] to identify which DEGs were significantly enriched in GO terms and metabolic pathways at Bonferroni-corrected P-value ≤ 0.05 compared with the whole-transcriptome background.

2.7. Quantitative Real-Time PCR (qRT-PCR) Validation.

Using the identical RNA/cDNAs for RNA-seq as templates, qRT-PCR was performed on a Bio-Rad CFX96 (Bio-Rad Laboratories) to verify the authenticity of transcriptomic profile expression patterns. The 10 μL reaction system contained TB Green®Premix Ex Taq™ II 5 μL, 10 μM primers each 0.2 μL, cDNA template 1 μL, and ddH₂O 3.6 μL. The amplification procedure was initially 95 °C for 30 s, followed by 40 cycles of 95 °C for 10 s and 60 °C for 30 s (two-step thermal cycling). Twenty-one candidate DEGs involved in various processes were randomly picked up as target genes, and the housekeeping gene TaActin was used as an internal control. The primers used for qRT-PCR are listed in Table S1. The relative level of expression was calculated using the formula 2^-ΔΔCT [36] (Livak and Schmittgen, 2001).

2.8. Weighted Gene Co-Expression Network Analysis (WGCNA)

Co-expression networks were constructed based on pairwise correlations of gene expression across all samples to find clusters (modules) of highly correlated genes and explore the connection of the modules with physiological traits and the hub genes [37]. After filtering out the genes with a low expression (FPKM < 2; CV > 0.15), 9924 genes were retained for the WGCNA. We introduced the appropriate soft threshold power β = 7 in this study according to Scale independence and Mean Connectivity curves produced by the software. Then, the adjacency was transformed into a topological overlap matrix (TOM), and average linkage hierarchical clustering was conducted according to the TOM-based dissimilarity measure, merged genes with a minimum size of 30, and Merge Cut Height 0.25 as a Module. Modules were defined as clusters of highly interconnected genes, and genes within the same cluster have high correlation coefficients among them. The total connectivity and intramodular connectivity were calculated with weighted and correlation coefficients in each module. Biologically meaningful co-expression modules were subsequently analyzed through the heatmap of the correlation coefficients between trait and module. Co-expression network analysis uses the top 30 gene nodes with connectivity, in which each node represents a gene and the connecting lines (edges) between genes represent co-expression correlations. Interaction network visualization for each module was performed using Cytoscape version 3.9.1 [38]. Hub genes have the highest Kme (Eigengene connectivity) value, indicating the most connections in the network.

2.9. Statistical Analysis

Statistical analysis was carried out using SPSS version 19.0 for Analysis of variance (ANOVA), followed by Duncan’s multiple range test (DMRT) to compare means among treatments if the ANOVA result was significant (p ≤ 0.05). Before ANOVA, the data were tested for normality by Chi-square analysis and transformed where necessary.

Since RRL is a ratio of means, while the measured net root growth values (x) were not paired between control and stress groups, the mean value and standard derivation of RRG were calculated according to Formulae:

$$\overline{RRL}=\frac{{\overline{x}}_{AlorMn}}{{\overline{x}}_{CK}};{SD}_{RRL}=\overline{RRL}\times \sqrt{{\left(\frac{{SD}_{AlorMn}}{{\overline{x}}_{AlorMn}}\right)}^2+{\left(\frac{{SD}_{CK}}{{\overline{x}}_{CK}}\right)}^2}$$

3. Results

3.1. Tolerance of Wheat Varieties to Al and Mn toxicity

Inhibition of root growth, as indicated by relative root length (RRL), is the initial and primary symptom of early Al toxicity. Aluminum tolerance was precisely the breeding target for creating NILs ES8 and ET8. Not surprisingly, the RRLs of aluminum-tolerant parent Carazinho and near-isogenic line ET8 were significantly higher than those of Egret and ES8. Interestingly, the change order of root length (Carazinho > ET8 > Egret > ES8) under Mn stress was the same as that under Al stress, and the change amplitude was similar (Figure 1).

Al and Mn toxicity decreased the chlorophyll level, as indicated by SPAD value, and the decreasing amplitude in Carazinho and ET8 were lower than those of Egret and ES8. There was no significant difference in the inhibitory effect of Mn and Al in Carazinho and ET8, while Mn showed a more substantial effect than Al in Egret and ES8 varieties. Furthermore, Mn toxicity symptoms like brown spots on the old leaves or crinkles on the young leaves were found only on leaves of Egret and ES8 after 4 d of Mn stress (Figure 2).

Collectively, nearly equal damage to wheat was caused by 25 µM Al and 267 µM Mn. The concentrations of active Al and Mn stresses screened from the preliminary experiments were suitable for further comparison of phenotypes and mechanisms of the two stresses. The consistency of tolerance to Al and Mn in wheat near-isogenic lines (ET8 and ES8) and their parents (Carazinho and Egret) confirmed the possibility of co-adaptation and breeding multi-resistant varieties.

3.2. Al and Mn Subcellular Distribution in Root Cell and Contents in Plants

After influx from the rhizosphere and entry into root cells, Al was mainly fixed in the cell wall, with a proportion variation of 82% to 93% among genotypes, and only a small part was distributed in cell saps. The proportions of Al fixed by cell wall in tolerant varieties were lower than those of sensitive varieties, possibly resulting in a low degree of damage of cell wall rigidification. By contrast, most of Mn in the root cell was distributed in the cell sap rather than the cell wall. There was a slight genotypic variation but no correlation to Mn tolerance (Figure 3).

Most of the Al absorbed by the wheat accumulated in the root, only a small part of Al (< 20%) was transported to the shoot, while Mn was transported to the shoot, and only a small amount of Mn (< 10%) stayed in the root (Figure 3). The root-to-shoot ratio (R/T) of Al and Mn contents did not differ significantly among cultivars (data not shown). Of note, due to the small biomass, the Mn and Al contents in the root, shoot, and whole plant of ES8 and ET8 were significantly lower than those of their parents Carazinho and Egret. There was no significant difference in concentration based on dry weight among cultivars (data not shown).

3.3. Al and Mn Effectiveness on Four Metal Ions Contents

The variance analysis results revealed that for every element (Ca, Mg, Fe, and Zn), the contents in wheat were significantly different among stress treatments (F = 27.336~83.02, F_0.01(2, 24) = 5.61), varieties (F = 10.9~116.62, F_0.01(3, 24) = 4.72), and their combinations (F = 19.69~52.90, F_0.01(11, 24) = 3.09). Ca, Mg, and Zn levels in wheat were dramatically decreased by both Al and Mn stress. The decreasing effects of Al and Mn toxicity on Ca content were roughly equivalent. However, the decreasing effects of Al stress on Mg and Zn contents were more substantial than those of Mn toxicity. Al toxicity significantly increased Fe contents, while Mn toxicity had the opposite effect. The significant differences in Ca, Mg, Fe, and Zn contents among varieties were not correlated to Al or Mn tolerance. The significance of multiple differences among stress combinations with varieties was labeled in Figure 4.

3.4. Wheat Root Transcriptome Profiling in Response to Al and Mn

Whole genome transcriptome sequencing analysis was applied to investigate the molecular responses of wheat roots to Al and Mn toxicity. Thirty-six cDNA libraries were constructed, combining three stresses, four varieties, and three replicates. These libraries yielded 1.8 G clean reads. The Q30 values of all samples exceeded 93.5% with a confidence of 99.9%, and the minimum mapping rate was over 90.55%, indicating that the sequencing was of high quality and qualification for the subsequent gene expression analysis.

A total of 115534 expressed genes (91609 known and 23925 new) and 193406 transcripts (113693 known and 79713 new) were assembled. Functions of 91313 genes and 113446 transcripts were finally annotated referencing KEGG, Swiss-Prot, Pfam, GO, COG, and NR libraries.

To validate the accuracy of gene expression of RNA-seq, qRT-PCR was performed using 21 randomly selected genes (Table S1) and identical cDNA templates. The results of qRT-PCR showed a very significant positive correlation with the transcriptome data, indicating that the transcriptome data were reliable (Figure 5).

3.5. Identification of Differentially Expressed Genes in Response to Al and Mn Toxicity

The FPKMs of genes expressed under Al and Mn toxicity were compared with control using the criteria|log2FC| ≥ 1 and P < 0.001. The DEGs from various varieties were pooled within stress treatments. A total of 5461 DEGs were identified under Al stress, comprising 3227 upregulated and 2234 downregulated DEGs (Figure 6A). Similarly, 4165 DEGs were identified under Mn stress, comprising 2617 upregulated and 1548 downregulated genes (Figure 6B).

The Venn diagram categorizes the 7699 genes responded by Al and Mn into eight sets, classified into four types based on expression patterns (Figure 6C, 6D, Table S2). Among these, 1,950 genes expressed in the same pattern under Al and Mn toxicity (1283 genes upregulated and 667 genes downregulated by both Al and Mn; Type I) were identified as possible candidate genes contributing to common toxicity or tolerant mechanisms. In contrast, seven genes differentially expressed under Al and Mn toxicity but in the opposite regulatory direction (Type Il), among which one gene annotated as solute carrier family 40 member I was downregulated by Al but upregulated by Mn, and six genes were upregulated by Al but downregulated by Mn, including three genes encoding bidirectional sugar transporter, one gene encoding ricin B-like lectin, one gene encoding cytochrome P450 and one gene encoding cysteine synthase. In addition, 3504 genes (1938 up, 1566 down) were exclusively regulated by Al (Type Ill), and 1950 genes (1333 up, 875 down) were exclusively regulated by Mn (Type IV). The latter three types of genes were studied to analyze the potential specific toxicity or tolerant mechanisms (Figure 6C 6D).

3.6. GO Functional Annotations and KEGG Pathway Analysis

The GO annotation results showed that binding (GO: 0005488) and catalytic activity (GO: 0003824) in the molecular function (MF) categories, membrane part (GO: 0044425), cell part (GO: 0044464), and organelle (GO: 0043226) in the cellular component (CC) categories; metabolic process (GO: 0008152), cellular process (GO: 0009987), and biological regulation (GO: 0065007) in the biological process (BP) categories were the most abundant DEGs of GO annotation terms. In transporter activity (GO: 0005215) term, 193 genes (88 up, 105 down) responded exclusively to Al, 161 genes (95 up, 66 down) exclusively responded to Mn, and 94 genes (49 up, 45 down) responded to both Al and Mn expressed in different trends. In transcription regulator activity (GO: 0140110) term, 145 genes (145 up, 35 down) responded exclusively to Al, 161 genes (109 up, 22 down) responded exclusively to Mn, and 94 genes (115 up, 14 down) responded to both Al and Mn, indicating that most transcription factors were upregulated to respond Al and Mn. In antioxidant activity (GO: 0016209) term, 90 genes (32 up, 58 down) exclusively responded to Al, 59 genes (15 up, 44 down) exclusively responded to Mn, and 94 genes (5 up, 38 down) responded to both Al and Mn (Figure 7).

The top 20 KEGG enrichment pathways with the smallest FDR values are displayed in Figure 8. The results showed that the genes upregulated by Al or Mn were enriched into five categories, including metabolism (M,14 pathways), genetic information processing (GIP, 4 pathways), environmental information processing (EIP, 2 pathways), organismal systems (OS, 1 pathway), and cellular process (CP, 1 pathway). In contrast, the genes downregulated by Al or Mn were mainly grouped into the metabolism (M) category, which contains 22 pathways. Genes upregulated by both Al and Mn were enriched in protein processing in endoplasmic reticulum (GIP,137 genes), spliceosome (32 GIP), plant-pathogen interaction (OS,31 genes), endocytosis (CP, 16 genes), and phenylalanine metabolism (M, 15genes). Genes downregulated by both Al and Mn were enriched in seven pathways, such as phenylpropanoid biosynthesis (M, 46 genes), starch and sucrose metabolism (M, 10 genes), and base excision repair (GIPs, 5 genes). Genes upregulated by Al were enriched in diterpenoid biosynthesis (M, 16 genes) and aminoacyl-tRNA biosynthesis (GIP,2 genes). Al downregulated genes were enriched into 12 pathways, such as fatty acid degradation (M, 10 genes), terpenoid backbone biosynthesis (M, 16 genes), and zeatin biosynthesis (M,10 genes). Genes upregulated by Mn were enriched into plant hormone signal transduction (EIP, 28 genes), alpha-linolenic acid metabolism (M, 8 genes), phenylalanine, tyrosine and tryptophan metabolism (M, 7 genes), and ABC transporter (6 genes) pathways. Genes downregulated by Mn were enriched into cutin, suberine and wax biosynthesis (M, 14 genes), benzoxazinoid biosynthesis (M,6 genes), and glycerolipid metabolism (M, 9 genes). Genes upregulated by Al and Mn were enriched in MAPK signaling pathway (EIP, 23 genes upregulated by Al and 50 genes upregulated by Mn) and brassinosteroid biosynthesis (M, 5 genes upregulated by Al and 4 genes upregulated by Mn), and genes downregulated by Al and Mn were enriched in galactose metabolism (M, 9 genes downregulated by Al and 9 genes downregulated by Mn). Genes were enriched in fatty acid (M, 8 genes upregulated by Al and 10 genes downregulated by Mn), and cysteine and methionine metabolism (M, 17 genes downregulated by Al and 23 genes upregulated by Mn) were oppositely regulated by Al and Mn.

3.7. Correlations between Physiological Traits and Expressed Module Eigengenes

Correlation analysis showed that RRL and SPAD were significantly positively correlated under Mn stress but not under Al stress. The Al contents in root, stem, and total plant were significantly positively correlated, while they are not correlated with the distribution ratio between cell wall and cell sap. Similarly, the Mn contents in root, stem, and total plant were also significantly positively correlated. Besides that, they all were negatively regulated with distribution in the cell wall and thus positively correlated with distribution in the cell sap. There was no significant correlation between Ca, Mg, Fe, and Zn variation amplitude under Mn stress. However, the variation of Mg was positively correlated with Ca and negatively correlated with Fe and Zn under Al stress (Figure 9).

A total of 9924 filtered genes in 36 samples were clustered and divided into 18 module eigengenes (ME) decorated with diacritical colors, respectively, based on a scale-free topological model β = 7 (Figure 10, Figure S1). The correlation coefficient and correlation between the module and trait were calculated to explore the genes and expression patterns involved in the physiological trait and tolerance and to find specific expression genes in samples. Five among eighteen MEs were notably correlated with at least one trait. MEmidnightblue (95 genes) was positively correlated with Mn traits (contents and distribution). MEblack (409 genes) was positively correlated with Al traits (contents and distribution). MElightcyan (93 genes) was positively correlated with Al traits and Fe content while negatively correlated with Mn traits. MEred (642 genes) was negatively correlated to RRL. MEturquoise (2075 genes) was positively correlated with RRL while negatively correlated with Al traits (Figure 10). The expression and annotation of genes in the above five modules are listed in Tables S3-1 to S3-5. The module is a highly interconnected set of genes that co-express to function and involve a biological process. However, the genes within a module are not necessarily differentially expressed genes. For instance, 2075 genes were clustered into MEturquoise; none were upregulated, and only a small portion was weakly downregulated by Al and Mn. RRL is the primary indicator for Al and Mn toxicity and tolerance screening. Forty-two genes in MEturquoise were directly related to root development such as putative auxin transporter-like protein (LAX3) related to lateral root formation, bifunctional L-3-cyanoalanine synthase (PALP) related to root hair cell development, germin-like protein (GLP) regulation of root development, protein root hair defective 3 (RHD3), root phototropism protein 3 (NPH3), and root-specific lectin (ChtBD) (Table S3-5). These genes did not respond to Al and Mn but provided additive information not unfolded by DEGs.

Connectivity analysis was conducted with the top 30 genes with high connectivity in each ME (MEmidnightblue, MEblack, MElightcyan, MEred, and MEturquoise) by Cytoscape to visualize the gene interaction regulation network (Figure 11).

The top five genes with a co-expression weight greater than 0.15, more connectivity, and higher Kme value in the interaction networks were screened. The candidate hub genes in the module were designated based on the annotated gene function and the correlation between ME and traits. The connectivity, expression, and annotation of candidate hub genes are shown in Table 1.

In MEmidnightblue, all the top five genes encoding nicotianamine synthase were upregulated exclusively by Mn. In MEblack, two genes function as zinc finger transcription factor, and three genes involved in intracellular trafficking and secretion were upregulated strongly by Al while marginally by Mn. In MElightcyan, two genes functioned as inorganic ion transporters and were involved in glutathione metabolism; one gene encoding TCP family transcription factor and one gene involved in protein processing in endoplasmic reticulum presented more connectivities. However, their expression was constitutive, not strongly regulated by either Al or Mn. This result provides new insights into the mechanisms underlying Al and Mn toxicity and tolerance beyond differentially expressed genes (DEGs). In MEred, three genes involved in genetic information processing and two genes encoded HSP20 were significantly upregulated by Al and Mn. In MEturquoise, two genes involved in proton and auxin transport and two genes involved in cell wall biogenesis and phenylpropanoid biosynthesis were slightly downregulated by Al and Mn. These genes may be a hub in their module and play a vital role in regulating the corresponding physiological trait.

4. Discussion

4.1. Wall Cell Biogenesis and Macromolecule Metabolism Respond Exclusively to Al Stress

The root tip is the first site for aluminum absorption, and aluminum toxicity inhibits the expansion and elongation of root cells and subsequent cell division [2]. Cell walls are the first barrier to Al uptake at the root apex. More than 80% of the total accumulated Al in plant roots is tightly bound to cell walls, with only a tiny fraction entering the cell sap (cytoplasm), while Mn is mainly distributed in cell sap (Figure 3). The primary cell wall components, pectin and hemicellulose, can bind Al; the binding of Al to cell walls decreases extensibility and cell elongation and increases rigidity [11]. Al stress significantly increases the contents of cell wall polysaccharides (pectin, hemicellulose I, and hemicellulose II) and the activities of pectin methylesterase and pectin demethylation at root tips, increasing the fixation of Al by cell walls [39]. The pectin methylation degree is determined by pectin methylesterase (PME) activity [15]. Xyloglucan is the most active Hemicellulose component in the cell wall; it is biosynthesized in the Golgi apparatus by a series of glycan synthases and glycosyltransferases (XTHs) before export to the wall [40]. XTHs are considered important factors in controlling the strength and ductility of the cell wall and are extremely sensitive to Al [41]. Becnel et al. (2006) found that at least twenty-one XTH members had strong expressions in the root, especially XTH14, XTH15, and XTH31 [42]. Yokoyama and Nishitani (2011) reported root-specific expression in XTH5, XTH12, XTH13, XTH14, XTH17, XTH18, XTH19, XTH20, XTH26, and XTH31 were differentially expressed in response to Al [43].

This study found that genes regulating cell wall modification were specifically induced by Al stress but not significantly by Mn. Four of the five genes encoding PME were upregulated, while one was downregulated, which may explain the increase in cell wall Al binding (Figure 12). Similarly, XTHs were differentially expressed with a divergent pattern: XTH22, XTH23, XTH25, XTH27, and XTH28 were upregulated, while XTH8, XTH16, XTH24, XTH26, XTH30, and XTH31 were downregulated. In addition, the NAC transcription factor, which directly binds to the XTH31 promoter region in Arabidopsis [44], was also highly upregulated in our study. In vivo, localization of XTH activity showed that Al greatly inhibited this enzyme activity within 30 min of exposure, which was concomitant with Al-induced callose deposition in roots [45]. Under Al stress, plants may increase resistance by reducing the production of xyloglucans through the downregulation of the XTH genes wherein XTH31 plays a significant role; the expression of other XTH genes may be involved in synergistic XTH31 expression or provide energy [41,45].

WGCNA results showed that MEturquoise and MEblack were strongly correlated with Al while weakly correlated with Mn. Numerous genes in the two MEs were involved in synthesizing and transporting cell wall components (Table S3-2, Table S3-5). For instance, peroxidase (POD) P7, encoded by the hub gene TraesCS7D02G369700, is mainly involved in clearing root cell ROS and synthesizing lignin. Another hub gene, TraesCS3B02G280000, annotated as V-type proton ATPase bundle B2, regulates the synthesis and transport of lignin monomers (Table 1) [46]. Lignin synthesis is a typical defense mechanism in plant stress response. Lignin accumulates in the cell wall of the root endothelial layer, inhibiting the entry of heavy metals into the xylem or their transport from vascular bundles to the aboveground. OsSTAR1 and OsSTAR2 OsSTAR1 and OsSTAR2, encoding a nucleotide-binding domain and a transmembrane domain of an ABC transporter, improve Al tolerance in rice through cell wall modification [10,47,48]. TraesCS3A02G129000 gene encoding ABC transporter C family member 3-like was exactly the hub gene in MEblack (Table 1).

The hardening of cell walls has a dual effect: On one hand, it inhibits the elongation of roots, and on the other hand, it prevents harmful metals from entering the cytoplasm. The modification of root cell wall properties can be reversible, so it is worth discussing whether cell wall stiffening is the cause or consequence of root growth inhibition.

4.2. Nicotianamine Synthesis Responds Exclusively to Mn Stress

Nicotianamine (NA) is a low-molecular-weight metal-bound ligand that maintains the homeostasis of Fe, Zn, and Mn. It protects cells from metal-induced damage, such as oxidative stress [49]. L-methionine is the precursor for NA biosynthesis, which binds with ATP to synthesize S-adenosylmethionine (SAM). Subsequently, three S-adenosylmethionine molecules are condensated to NA by nicotianamine synthase (NAS). NA is finally converted into 2’-deoxymugineic acids (DMA) and other mugineic acids (MAs) under the catalysis of nicotianamine aminotransferase (NAAT) [50]. Overexpression of MxNAS2 and MxNAS3 genes from M. xiaojinensis increases the concentration of NA-Mn in flowers and leaves of tobacco and Arabidopsis, respectively [51,52]. Yellow strip-like protein (YSL) is involved in the uptake, transport, and relocation of metal complexes such as Mn-NA, Mn-HMA, and Mn-DMA to maintain Mn ion balance in plants [53,54].

The results included nine genes encoding NAS, four genes encoding NAAT, and two genes encoding YSL, specifically and strongly upregulated by Mn stress but not significantly regulated by Al (Figure 13). For instance, the expression of NAS4, NAS1, NAS2, and NAATA increased 3040, 420, 342, and 577-fold under Mn stress. In addition, the YSL2 gene (TraesCS6D02G223000), which transports Mn-NA and Mn-DMA [54,55,56], was also upregulated exclusively by Mn (Figure 14).

WGCNA showed that the genes encoding NAS, NAAT, and YSL2 were all significantly upregulated by Mn in MEmidnightblue, along with three genes encoding heavy metal-associated domain (HMA) and one gene encoding ZIP were identified (Table S3-1). Most importantly, five genes encoding NAS were all included in the top 30 genes with the highest connectivity (Table 1), playing a hub role in the regulatory network to Mn stress.

Increasing the synthesis of NA and its derivatives and enhancing the long-distance transport of Mn-MAS is a specific strategy to tolerate Mn.

4.3. Metal Ions Transportation and Accumulation under Al and Mn Stress

The Al absorption site is located at the root tip, and Al induces a series of transporter genes to coordinate the absorption, transport, and redistribution of Al. NRAMP aluminum transporter 1 (NRAT1) plays a vital role in rice Al tolerance by reducing the level of toxic Al in the root cell wall and transporting Al ions (not Al-citrate complex) into the root cell, where they are ultimately sequestered in the vacuole [28,57]. The ALS protein family is responsible for transporting excess Al to the vacuole for isolation and thus enhancing tolerance to Al in rice [58], tea [59], and buckwheat [60]. In this study, Under Al stress, two nodulin 26-like intrinsic protein (NIP) genes TaNIP1;1 (TraesCS7B02G122600) and TaNIP1;2 (TraesCS7A02G215700) were significantly upregulated by Al stress. It may promote the transport of malate-Al complex (not Al ion) from the cell wall to the symplasm and further transport to the leaf vacuole for segregation [14,47,61]. WGCNA showed that MEblack and MElightcyan were significantly correlated with the content and distribution of Al. Some genes upregulated in MEblack may be involved in Al absorption, including those annotated as ABC transporter and one gene encoding NRAMP (Table S3-2).

Manganese is essential for plant growth, with a narrow concentration window between deficiency and toxicity. Therefore, Mn transport systems should be tightly regulated in plants. Plants must deliberately control the uptake, transportation, and distribution of manganese to prevent toxicity caused by excessive amounts of this element. Transporters play a vital role in these processes and thus help plants to tolerate Mn toxicity. Mn absorption occurs in the mature root zone and can be easily transported through the xylem to the shoot [62]. The Arabidopsis root plasma membrane-localized transporter natural resistance-associated macrophage protein (AtNRAMP1) mediates Mn²⁺ absorption. Two members of the zinc, iron-regulated transporter protein (ZIP) family, AtZIP1 and AtZIP2, are involved in the transport of Mn from root to stem [63,64]; OsYSL2 is responsible for long-distance transport and distribution of Mn [56]. In this study, the genes of TaYSL2 (TraesCS6D02G223000), TaZIP1 (TraesCS2A02G424200), and TaZIP2 (TraesCS6A02G158700) were significantly upregulated under Mn stress, and thus increased Mn absorption (Figure 14). Mn shares some ion channel proteins with other metal ions. For instance, OsNRAMP5, a plasma membrane protein, is involved in constitutive Fe and Mn uptake and inducive transport during flowering and seed development [65]. In the present study, the expression of reported multifunctional transporter genes, such as CAX2, CCX3, MTP9, MTP10, MTP11, MTP11, NRAMP3, NRAMP5, and YSL6 [63], were not regulated by Mn. The diversity of genes involved in Mn uptake and transport and the versatility of genes lead to the complexity of the regulation of Mn ion balance.

Al and Mn inhibit root growth, subsequently leading to the deprivation of essential nutrient elements. Due to its higher charge, the Al ion inhibits the absorption of bivalent cations (Ca²⁺, Mg²⁺, Fe²⁺, and Zn²⁺) more strongly than other toxic ions [66]. It removes bivalent cations from the plasma membrane and cell wall [67]. High concentrations of Mn²⁺ compete with other nutrient cations, thus hindering their uptake [63]. On the other hand, roots can improve tolerance to Al and Mn toxicity stress by controlling the uptake of mineral nutrients and maintaining ion homeostasis, including calcium [68,69], magnesium [70,71], zinc, and iron [63,72,73].

Our results showed that Al and Mn stress decreased Ca, Mg, and Zn contents, the absorption and accumulation of Fe by Al, while Mn exerted the opposite effect. Al and Mn significantly downregulated one gene (TraesCS5D02G030700) encoding the cation/calcium exchanger (CCX1), contributing to the decrease in Ca content. Plants utilize calmodulin-activated Ca²⁺ pumps (ACA) at the plasma membrane, endoplasmic reticulum, and vacuole [74]. The ACA7 (TraesCSU02G003600) and ACA5 (TraesCS5D02G284100) might play an essential role in altering the production of cellular Ca²⁺ signals and thus mediate adaptive responses (Figure 14).

Al and Mn shared hydrating ionic radii with Mg, competing for binding sites on apoplasm and transporters in the plasma membrane [75]. Aquaporin (PIP) regulates water and Mg²⁺ transport [76]. Two PIPs (TraesCS5A02G100000 and TraesCS7B02G002000) were significantly downregulated, related to a decrease in Mg²⁺ absorption. The decrease of Mg content under Al and Mn stress may be the reason for the decrease in chlorophyll content (SPAD). Overexpression of MGT1 can enhance Al tolerance [71,77]. The expression of one MGT gene (TraesCS3A02G380600) was upregulated under Al stress but not Mn due to an adaptive response related to the regulation of membrane potential balance by H⁺-ATPase activity under Al stress [77] (Figure 14).

Fe is an essential trace element for plant growth and development. However, free Fe in cells produces toxic ROS through the Fenton reaction [78], necessitating strict control of iron homeostasis through various strategies. Al stress increased Fe absorption in tea plants by upregulating Fe ion transporter genes, including FRO, YSL, and IRT, maintaining Fe homeostasis [59]. Mn and Fe have antagonistic solid effects, and excessive Mn²⁺ takes advantage of Fe on the binding of Fe/Mn co-transporters [63,79]. Most of the CDF/MTP and VIT families can transport Fe²⁺ and Mn²⁺, and excess Mn²⁺ blocks the absorption of Fe²⁺ [80,81]. In this study, three VIT family genes and one YSL5 gene were significantly downregulated under Al and Mn stresses, reducing Fe absorption. One IRT1 gene (TraesCS4D02G017600) and one MGP gene (TraesCS3D02G257900) were upregulated by Al stress while downregulated by Mn stress (Figure 14). These results indicate that Al stress increased Fe²⁺ absorption by upregulating IRT and MRP. MTP4 (TraesCS4A02G077100) and MTP7 (TraesCS3D02G038000) were downregulated and upregulated by Mn stress, respectively, which affected the absorption and transport of Mn and Fe.

High active Al³⁺ competes with Zn²⁺ for binding sites on the root plasma membrane, and excess Mn²⁺ shares a similar transporter with Zn²⁺, resulting in plant inhibition of zinc absorption [73,82]. Heavy-metal-associated domain (HMA) and ZIP proteins are involved in the absorption and transport of Zn [83]. Our study indicated that seven genes from the HMA family were significantly downregulated under Al and Mn stress. The expression of genes from the ZIP family was diverse, with the significant upregulation of ZIP1 (TraesCS2A02G424200) and the downregulation of ZIP8 (TraesCS2B02G533800) under Al and Mn stress. In addition, ZIP2 (TraesCS6A02G158700) was upregulated by Mn stress but downregulated by Al, which is favorable for the absorption of Mn and Zn (Figure 14).

In MEturquoise (Table S3-5), the genes encoding the calcium load-activated calcium channel (CAEEL), Fe transporters VIT and YSL, and Zn transporter HMA were downregulated by Al and Mn. In contrast, Zn transporter genes (ZIP) in the MEred were upregulated (Table S3-4). MEturquoise was positively correlated with RRL, while MEred was negatively correlated with RRL. Reducing relative root length may decrease absorption area and gene expression levels, ultimately inhibiting nutrient uptake.

4.4. Phenylpropanoid Biosynthesis under Al and Mn Stresses

The KEGG results indicated that genes involved in phenylalanine metabolism commonly responded to Al and Mn, and most genes involved in phenylpropanoid metabolism were downregulated by both Al and Mn; some genes involved in flavonoid metabolism specifically responded to Al (Figure 8). There is the same expression trend of genes involved in the metabolism pathway from phenylalanine to sinapoy-choline and scopolin under Al and Mn, with the downregulation of seven genes encoding phenylalanine ammonia-lyase (PAL), three genes encoding putrescine hydroxycinnamoyltransferase (PHT), one gene encoding caffeoylshikimate esterase (CSE), two genes encoding tricetin 3',4',5'-O-trimethyltransferase (COMT), three genes encoding Cytochrome P450 (CYP), and three genes encoding Anthocyanin 3'-O-beta-glucosyltransferase (TOGT), and upregulation of one gene (TraesCS3A02G329200) encoding serine carboxypeptidase-like (SCPL), except for one gene (TraesCS7A02G398900) encoding CYP45073A and one gene (TraesCS5B02G522400) encoding TOGT1 were exclusively upregulated by Al. Two genes (TraesCS6D02G136600 and TraesCS6A02G147300) encoding cinnamoyl-CoA reductase (CCR) were upregulated by Al and Mn. Among seven genes encoding beta-glucosidase (BGLG), three were downregulated by Al and Mn, BGLG31 (TraesCS5D02G302700) and BGLG6 (TraesCS4A02G056400) showed significant upregulation by Al, BGLB12 (TraesCS2A02G329000) and BGLB3 (TraesCS5D02G404700) were upregulated by Mn but downregulated by Al. One gene (TraesCS3D02G388100) encoding 4-coumarate-CoA ligase (4CL) and one gene (TraesCS7D02G328300) encoding cinnamyl alcohol dehydrogenase (CAD) were exclusively up and downregulated by Al, respectively (Figure 14).

Phenylpropanoid biosynthesis is an important pathway to produce plant secondary metabolites, such as phenylpropanoids, flavonoids, and diterpenoids, which are secreted in response to metal ions [84,85,86]. PAL catalyzes phenylalanine cinnamic acid, the first rate-limiting enzyme in the phenylalanine biosynthesis pathway [87]. Our results showed that genes involved in the pathway from cinnamic acid to sinapoy-choline and scopoline were downregulated, while PAL and CCR were upregulated by both Al and Mn. Similar to our study, genes encoding PAL and CCR were upregulated in phenylpropanoid metabolism and are activated to counteract Mn stress in Stylosanthes [88] and Arabis paniculate [89]. In Populus euphratica, the content of phenylalanine and the antioxidant enzyme activity increases under Mn stress [90]. Phenylpropanoid metabolism exhibits different expression patterns in the roots and leaves under Mn stress in Stylosanthes. The expression of PAL1 and chalcone-flavanone isomerase (CFI) was upregulated in the leaf, while PAL2 and chalcone synthase (CHS) were downregulated in roots [91]. It was also found that phenylpropanoid biosynthesis and phenylalanine metabolism played a critical role in defense against Aluminum [92]. Lettuce protected roots from oxidative damage by upregulating phenylalanine ammonia-lyase (PAL), cinnamate 4-hydroxylase, and 4-coumarate coenzyme A ligase [93]. These results suggest that phenylpropanoid biosynthesis enhanced antioxidant capacity for tolerance of both Al and Mn.

In addition, Al specifically upregulated 4CL and CYP, which are involved in the synthetic precursors of flavonoid metabolism, p-coumaroyl-CoA and cinnamoyl-CoA [84]. The KEGG results showed Al exclusively downed some genes involved in flavonoid metabolism, and two genes (TraesCS2D02G582000 and TraesCS4A02G003100) encoding tryptamine benzoyltransferase 1 (TBT1) and one gene (TraesCS2A02G482200) encoding Anthocyanidin reductase (ANR) were exclusively upregulated (Table S2-3). The flavonoid biosynthesis site is located at the upper end of the surface cells in plants’ root growth area and participates in internal and external Al detoxification by forming solid complexes with toxic Al ions [94]. Alfalfa increased resistance to Al stress by increasing flavonoid accumulation in roots and secretion from root tips [95]. Our results suggest wheat tolerance of Al by promoting flavonoids to chelate Al ions.

Figure 15. The expression levels of DEGs involved in phenylpropanoid metabolism under Al and Mn stresses. PAL: Phenylalanine ammonia-lyase; 4CL: 4-coumarate--CoA ligase; CCR: Cinnamoyl-CoA reductase; BGLG: Beta-glucosidase; CYP: Cytochrome P450; CAD: Cinnamyl alcohol dehydrogenase; TOGT: Anthocyanin 3'-O-beta-glucosyltransferase; PHT: Putrescine hydroxycinnamoyl transferase; CSE: Caffeoylshikimate esterase; COMT: Tricetin 3',4',5'-O-trimethyltransferase; SCPL: Serine carboxypeptidase-like protein.

Figure 15. The expression levels of DEGs involved in phenylpropanoid metabolism under Al and Mn stresses. PAL: Phenylalanine ammonia-lyase; 4CL: 4-coumarate--CoA ligase; CCR: Cinnamoyl-CoA reductase; BGLG: Beta-glucosidase; CYP: Cytochrome P450; CAD: Cinnamyl alcohol dehydrogenase; TOGT: Anthocyanin 3'-O-beta-glucosyltransferase; PHT: Putrescine hydroxycinnamoyl transferase; CSE: Caffeoylshikimate esterase; COMT: Tricetin 3',4',5'-O-trimethyltransferase; SCPL: Serine carboxypeptidase-like protein.

4.5. Signaling in Root under Al and Mn Stresses

According to the KEGG results, the majority of genes involved in signaling, including those related to the MAPK signaling pathway, plant-pathogen interaction, and plant hormones like IAA, ETH, and JA, exhibited a similar expression pattern. However, the specificity of ZT and BR in response to Al and Mn was noticeable.

Genes expressed commonly in plant hormones (Figure 8), including two genes (TraesCS2B02G210600 and TraesCS2D02G191800) encoding indole-3-acetic acid-amido synthetase GH3, one gene (TraesCS3D02G411300) encoding auxin-responsive factor (ARF), two genes (TraesCS2D02G391400 and TraesCS1A02G328800) encoding ethylene-responsive transcription factor (ERF), and six genes encoding jasmonate ZIM domain protein (JAZ) were upregulated. Two genes (TraesCS3A02G348400 and TraesCS3D02G342000) encoding abscisic acid receptor PYL were downregulated by both Al and Mn, with a pronounced downregulation by Al than Mn. While some precisely response to Al or Mn, three genes involved in brassinosteroid (BR) were exclusively upregulated by Al, while one gene (TraesCS5A02G214800) encoding CYP92A6 was exclusively upregulated by Mn. Furthermore, three genes involved in zeatin (ZT) were exclusively downregulated by Al, including two genes (TraesCS3A02G263300 and TraesCS3D02G263000) encoding adenylate isopentenyl transferase 1 (IPT), and one gene (TraesCS2B02G013000) encoding protein NRT1/ PTR family 8.5 (PTR) (Figure 16).

IAA and JA affected the antioxidant defense mechanisms. They altered the cellular redox homeostasis for Al and Mn, showing the commonality of plant hormones in response to Al and Mn stress [96]. In addition, plants respond to both stresses in some specific mechanism. The primary lesion of Al toxicity is the rapid inhibition of root elongation, plant hormones primarily reacting to Al stress in the root. Genes involved in BR and ZT were specifically reposed to Al, which might have contributed to root growth and organic acid. The activated ethylene signal by Al promoted IAA or cytokinin (CTK) accumulation and led to a modification of cell wall structure through auxin-responsive factors (ARFs). It was reported that BZ signaling is involved in the expression of ALMT1 and STOP. It also stimulated the synthesis of the jasmonic acid signaling to modulate microtubule (MT) polymerization to regulate root growth [15]. Auxin-regulated GmMATE and H⁺-ATPase enhance Al tolerance by modifying cell walls and promoting citrate exudation [97]. Abscisic acid (ABA) was found to mediate the Al-induced citrate exudation in soybean [98] and ethylene-mediated JA signal-regulated ALMT1-controlled malate exudation [99,100]. These results suggest that plants participating in Al tolerance mainly focus on the roots. While the changes in plant hormone content were influenced by Mn accumulation, moderate Mn increased IAA, ZT, and 6-BA (6-benzylaminopurine) contents but decreased ABA content in the leaves. High Mn accumulation significantly increased leaf JA and ABA contents, which inhibited plant growth and stimulated stress tolerance [101]. Mn enhanced resistance to Mn toxicity based on IAA and ABA synthesis by regulating growth and scavenging ROS for antioxidant and Mn transporter [102].

Most genes involved in the MAPK signaling pathway and plant-pathogen interaction expressed the same trend under Al and Mn stresses. Six genes encoding calmodulin-like proteins (CML) and one gene (TraesCS4B02G327800) encoding respiratory burst oxidase homolog protein (RBOH) were significantly upregulated by Al and Mn. Among three genes encoding mitogen-activated protein kinase (MAPK), M2K9 and M3K17 showed significant upregulation, while Al and Mn stresses downregulated M3K18. Moreover, one gene (TraesCS7D02G161200) encoding a pathogenesis-related protein (PRP) and MEKK1 was exclusively upregulated by Mn (Figure 16).

Signaling, such as Ca²⁺, ROS, and NO, is a late response to Al stress. Al promoted the binding of Ca²⁺ to CML24, leading to malate exudation by regulating the expression of STOP1, WRYK46, and ALMT1 [103]. The MAPK positively regulates STOP1, thus conferring Al resistance [104]. Similarly, the response of these signals is vital for Mn. Mn²⁺ increased catalase activity by regulating mitogen-activated protein kinases (MAPK) and calmodulins in the presence of calcium to reduce Mn toxicity [105]. The production of one of the ROS, hydrogen peroxide (H₂O₂), is often generated contemporaneously with NO, regulated by CML and RBOH, respectively. NO was shown to play a familiar role in reducing Al toxicity by regulating the ascorbic acid glutathione cycle and cell wall modification in wheat roots [39,106] and alleviating Mn toxicity by preventing oxidative stress in rice leaves [107].

4.6 Transcription Factors Regulate Wheat Tolerance to Al and Mn Stress

Go annotation analysis showed that transcription factors (TFs) were predominantly upregulated by Al and Mn (Figure 7). Transcription factors are the terminal points of stress signal transduction and molecular switches for downstream gene expression [108]. Numerous TFs, such as WRKY, MYB, GATA, bZIP, bHLH, and ERF, have been reported to regulate Al and Mn stress responses [25,109,110,111,112]. For instance, sensitive to proton toxicity (STOP1) and its rice homolog Al resistance transcription factor (ART1) are both members of the C2H2 zinc finger family, which enhances Al tolerance by regulating ALMT1 [113,114]. Overexpression of GsERF1 enhances aluminum tolerance in Arabidopsis thaliana through ethylene-mediated pathways and ABA signaling pathways [115]. Similarly, ZmbHLH105 improves tolerance to manganese stress by regulating antioxidant mechanism-mediated ROS clearance and Mn/Fe-related transporter expression in plants [109]. Some transcription factors that regulate genes associated with aluminum (Al) tolerance may also have a similar or opposite role in manganese (Mn) tolerance. For instance, C2H2 transcription factors activate Al-tolerant genes glutamate dehydrogenase 1 (GDH1), pectin methylesterase inhibitor (PMI), malic enzyme (ME), aluminum-activated malate transporter (ALMT), multidrug and toxic compound exclusion protein (MATE), and tonoplast dicarboxylate transporter (TDT) [116,117,118]. These homologous genes were upregulated by Mn toxicity [89]. WRKY46 is a negative regulator of ALMT1 in Arabidopsis Thaliana, and the destruction of WRKY46 leads to increased malate secretion, resulting in Al tolerance [112]. However, WRKY gene transcripts in stylo are enhanced by Mn stress [89], suggesting that there are conserved and divergent regulatory networks mediated by WRKY-mediated gene expression under Al and Mn toxicity.

Expression analysis showed that 640 TFs were significantly regulated by at least one stress, covering 29 TFs families (Table 2). The MYB family accounted for the most, followed by WRKY, AP2, bHLH, and NAC families. The differentially expressed TFs from BES1, GRF, HD-ZIP, RAV, TCP, and NF-YA families were exclusively regulated by Al. In contrast, two TFs from calmodulin-binding transcription activators (CAMTA) and ethylene-insensitive 3-like (EIL) families were exclusively regulated by Mn, indicating the specificity of these TFs in sensing stress and regulating downstream genes. A total of 455 differentially expressed TFs were found under Al stress, comprising 291 upregulated and 82 downregulated. Similarly, 373 TFs were differentially expressed under Mn stress, including 414 upregulated and 212 downregulated. The upregulated TFs were far more than the downregulated by either Al or Mn, suggesting that wheat tends to activate TFs to adapt to toxicity. Numerous differentially expressed TFs were found under only one stress, with 267 TFs exclusively regulated by Al and 185 by Mn. 188 differentially expressed TFs, including 156 upregulated and 32 TFs downregulated, were identified under two stressed in the same trend, demonstrating that most co-expressed TFs function similarly under Al and Mn stress.

5. Conclusions

The current study analyzed the physiology and transcriptomics of roots from wheat seedlings tolerant to Al and Mn stress. Results showed that Al was mainly fixed in the root cell walls, regulated by PEMs and XTHs. Mn was distributed primarily on the cell saps and then transferred to shoot by NA and YSL. Expression of metal cation transporter genes was altered by Al and Mn, leading to reduced absorption of Ca, Mg, and Zn. Additionally, Mn decreased Fe absorption, while Al had the opposite effect. The identified DEGs associated with commonality and specific pathways based on the KEGG enrichment analysis might represent an adaptive strategy of wheat to cope with Al and Mn. Numerous genes related to TFs, such as MYB, WRKY, and AP2 families, were activated to tolerate Al and Mn stress. In addition, WGCNA identified the hub genes PODP7, VATB2, and ABCC3 for Al tolerance and NAS for Mn tolerance. For these reasons, further study is recommended to pyramid genes to breed multi-tolerant varieties, understand the interaction of the TFs with other proteins, and regulate the target genes.