Aluminum Responsive Genes in Flax (Linum usitatissimum L.)

Flax (Linum usitatissimum L.) is a multipurpose crop which is used for the production of textile, oils, composite materials, pharmaceuticals, etc. Soil acidity results in a loss of seed and fiber production of flax, and aluminum toxicity is a major factor that depresses plant growth and development in acid conditions. In the present work, we evaluated gene expression alterations in four flax genotypes with diverse tolerance to aluminum exposure. Using RNA-Seq approach, we revealed genes that are differentially expressed under aluminum stress in resistant (Hermes, TMP1919) and sensitive (Lira, Orshanskiy) cultivars and selectively confirmed the identified alterations using qPCR. To search for differences in response to aluminum between resistant and sensitive genotypes, we developed the scoring that allowed us to suggest the involvement of MADS-box and NAC transcription factors regulating plant growth and development and enzymes participating in cell wall modifications in aluminum tolerance in flax. Using Gene Ontology (GO) enrichment analysis, we revealed that glutathione metabolism, oxidoreductase, and transmembrane transporter activities are the most affected by the studied stress in flax. Thus, we identified genes that are involved in aluminum response in resistant and sensitive genotypes and suggested genes that contribute to flax tolerance to the aluminum stress.


Introduction
Among abiotic stresses, aluminum (Al) toxicity is a major constraint for crop production in acid soils worldwide [1]. In acidic conditions, the mineral form of Al dissolves to release the soluble Al 3+ form, which is capable of crossing the plant membranes and is highly toxic to plants that even micro concentrations can inhibit root growth within minutes or hours in many agricultural plant species [2][3][4][5]. Al negatively affects cell elongation and division, uptake and transport of nutrients, and Ca 2+ homeostasis and disturbs the structure and function of the plasma membrane, cell wall, and chromatin [6][7][8][9][10].
The mechanisms of resistance to Al are diverse in plants and could be divided into exclusion, which decreases the amount of phytotoxic Al 3+ in the cells and internal tolerance, which reduces Al toxicity in root and shoot symplast [6,[11][12][13][14][15][16][17][18][19]. Excretion of Al detoxifying organic acids (OAs) to the apoplast or rhizosphere is the most common variant of exclusion, in which genes encoding the OA transporters, including aluminum-activated malate transporter 1 (ALMT1) and members of the multidrug and toxic compound extrusion (MATE) citrate transporter gene family, are involved [13,20]. Exudation of mucilage also binds Al ions and results in the exclusion of Al and its detoxification [21][22][23]. Mechanism of Al tolerance encompasses processes that result in chelation of cytosolic Al 3+ with organic ligands, sequestration into the vacuole, or transport of Al 3+ to less-sensitive regions of the plant [24][25][26]. Recovery from damages following exposure to Al toxicity is mediated by the detoxification of the reactive oxygen species (ROS) through ROS-detoxifying enzymes, such as glutathione S-transferases, peroxidases, and superoxide dismutases that confer Al tolerance [20,27,28]. Significant genetic diversity in Al resistance or tolerance was found in crops that enable development of improved cultivars, which maintain yield on acid soil [29,30].
For flax (Linum usitatissimum L.), a crop grown worldwide for fiber and seeds and also known for its health benefits [31,32], the mechanisms for Al tolerance are little known, although Al toxicity in acid soils is a serious problem for cultivation and rich harvest of flax [33]. Therefore, it is imperative to understand the molecular mechanisms underlying Al tolerance in flax. High-throughput sequencing methods are extensively used for evaluation of expression alterations of protein-coding genes or miRNAs in flax under diverse stresses [34][35][36][37][38][39][40][41]. We have previously shown that UDP-glycosyltransferases, glutathione S-transferases, and Ca 2+ /H + -exchanger CAX3 are involved in flax tolerance to aluminum [42,43]. miR319, miR390, and miR393 also participate in aluminum response via regulation of growth processes [44]. However, Al tolerance mechanisms are complex and could involve diverse strategies. In the present work, we conducted comprehensive gene expression analysis of flax cultivars grown under control and Al-treated conditions to reveal genes that participate in Al response and to identify differences in gene expression alterations between resistant and sensitive to aluminum genotypes.

Plant Material.
A laboratory method for evaluation of flax tolerance to aluminum, whose results have high correlation with field assessment, was taken as the basis for creating aluminum treatment conditions [45]. Resistant and sensitive cultivars for the study were chosen on the basis of our field and laboratory experiments, in which productivity indexes and root length were used for assessment of flax genotype tolerance to soil acidity and aluminum [33]. For the present study of gene expression alterations under shortterm (4 h) aluminum exposure (500 M AlCl 3 ) at low pH (4.5), resistant (Hermes and TMP1919) and sensitive (Lira and Orshanskiy) to aluminum stress flax cultivars were used. Flax seeds were germinated in Petri dishes on filter paper and then were transferred into 50 mL Falcon tubes with filter paper soaked in 0.5 mM CaCl 2 pH 4.5 for 24 h. Seedlings in control conditions were then grown for 24 h more, while seedlings in stress conditions were grown for 20 h followed by addition of 500 M AlCl 3 for 4 h. After that, root tips, 8-10 mm in length, were collected and immediately frozen in liquid nitrogen. The samples were stored at -70 ∘ b. RNA was extracted using RNA MiniPrep kit (Zymo Research, USA) and then used for transcriptome sequencing on HiSeq2500 (Illumina, USA) with paired-end 100-nucleotide reads in two biological replicates as described in our previous work [42].

Identification of Differentially Expressed
Genes. Transcriptome assembly and annotation were performed using Trinity and Trinotate as described earlier [42]. The reads of each cultivar (Hermes, TMP1919, Lira, and Orshanskiy) were mapped to the assembled transcripts of Hermes using bowtie2 [46] and quantified using RSEM [47]. Then, read count data were analyzed using edgeR package [48]. Genes with read counts per million (CPM) below 1.5 were filtered out. Normalization using the TMM method was performed and genes with expression alterations under aluminum treatment were identified within the following groups: (1) pool of all analyzed flax genotypes; (2) pool of resistant cultivars (Hermes and TMP1919); (3) pool of sensitive cultivars (Lira and Orshanskiy); (4) individual genotypes (Hermes, TMP1919, Lira, or Orshanskiy).
To evaluate expression alterations for each transcript, expression level fold change (FC) was calculated for each gene as the ratio of CPM under Al 3+ exposure to CPM under control conditions.
To identify genes with diverse expression alterations in resistant and sensitive to aluminum genotypes, we calculated delta-score ( Δ ) that takes into account (1) consistency of gene expression changes within each group (resistant or sensitive cultivars), (2) magnitude of these changes, and (3) differences between groups. To proceed this way, we derived consistency scores for resistant and sensitive cultivars, which correspond to the first component: where log is the binary logarithm of FC, and i = 1, 2 (number of resistant genotypes). Hence, if aluminuminduced gene expression changes are unidirectional and FC values are close to each other, res. would have the greatest value. In the similar way, we calculated the consistency score for sensitive cultivars ( sens. ). Finally, delta-scores ( Δ ) were calculated as Here, the 1 st multiplier reflects the consistency of gene expression changes, the 2 nd one-their maximum magnitude, and the 3 rd one-differences between groups. We introduced 2-fold greater weight for the 3 rd component because of its prime importance.

Gene Ontology Analysis.
The gene set enrichment analysis (GSEA) based on Gene Ontology (GO) data was performed using Goseq package (http://bioconductor.org/ packages/release/bioc/html/goseq.html). Analysis of top lists of up-and downregulated genes was carried out for the pool of all analyzed genotypes, the pool of resistant cultivars, the pool of sensitive cultivars, and individual genotypes for identification of enriched GO terms. Heatmaps illustrating gene expression alterations in selected GO terms were created using R package pheatmap. Weighted gene correlation network analysis was carried out using R package WGCNA [51] for genes with CPM > 20 in at least 6 samples.

The Effects of Aluminum Exposure at Low pH on Flax
Plants. Soil acidity and Al 3+ ions result in depression of flax plants, inhibition of growth, changes in height, fiber mass, and seed productivity both in resistant and sensitive cultivars. However, the degree of changes varied between genotypes [33]. For the present study, we used flax cultivars with contrast extent of phenotypic changes in acid soils with high aluminum content (pH 4.0, Al content -11.07 mg/100 g of soil). Results of field experiments for Hermes, fV`1919, Lira, and Orshanskiy cultivars under control conditions and low pH with Al treatment are presented in Table 2. Cultivars Hermes and fV`1919 showed greater tolerance compared to Lira and Orshanskiy that was especially noticeable for fiber mass and a number of seed pods.
In laboratory experiments, inhibition of root growth was observed in flax plants under aluminum exposure at low pH (pH 4.5, 500 M AlCl 3 ). Negative effects were more pronounced in sensitive cultivars, where 50-60% reduction in root length was revealed after 5 days, while in resistant cultivars reduction was less than 30% [45].

Genes with Differential Expression under Al 3+ Exposure.
As shown in our previous works, even after 4 h of Al 3+ exposure when phenotype changes were not yet noticeable, significant gene expression alterations occurred in flax plants [42,44]. These results are consistent with the studies on other plant species, for many of which short response time to aluminum was revealed [19]. Therefore, in the present work, we focused on the effects of short-term aluminum treatment and evaluated gene expression alterations after 4 h of aluminum exposure in resistant and sensitive to the studied stress flax genotypes based on transcriptome sequencing data. In this regard, in contrast to our previous study [42], we compared gene expression levels in flax plants under control conditions and after 4 h of Al 3+ exposure in the pool of all studied genotypes, groups of resistant and sensitive plants, or individual cultivars. Results of our analysis are presented in Supplementary Materials for the pool of all studied genotypes (S1 table), pool of resistant (Hermes and TMP1919, S2 Table) and pool of sensitive (Lira and Orshanskiy, S3 Table) genotypes, and individual cultivars (S4, S5, S6, and S7 Tables for Hermes, TMP1919, Lira, and Orshanskiy respectively), where the genes are listed in the order of decreasing statistical significance of expression alterations.
We performed the search for differences in transcriptomic response to aluminum stress between resistant and sensitive genotypes to identify genes potentially involved in aluminum tolerance. We developed the scoring that takes into account gene expression changes within resistant and sensitive groups of genotypes, their consistency, and differences between groups. This allowed us to find out transcripts with diverse expression alterations between groups of resistant and sensitive to aluminum flax genotypes but similar alterations within resistant or sensitive groups. Results of the analysis are presented in S8 Table. Transcripts were sorted by decreasing delta-score that reflects differences in Al 3+ -induced transcriptomic response between resistant and    sensitive groups of genotypes. Within upregulated transcripts in resistant genotypes compared to sensitive ones, transcripts encoding Agamous-like MADS-box protein AGL62, polygalacturonase, NAC domain-containing protein 100, protein OSB1, GDSL esterase/lipase, beta-glucosidase 11, (+)neomenthol dehydrogenase, osmotin-like protein OSML13, and 1-aminocyclopropane-1-carboxylate oxidase homolog 1 were in the top.
First of all, our attention was attracted by genes that were upregulated only in resistant cultivars under Al 3+ exposure. In the top of S8 Table, two TFs are located: Agamous-like MADS-box protein AGL62 and NAC domain-containing protein 100. MADS-box genes control numerous aspects of plant development and are involved in stress response [52][53][54][55]. NAC TFs also participate in stress response and their overexpression in transgenic plants enhances tolerance to abiotic stresses and promotes development of lateral roots [56][57][58][59][60][61]. Role of NAC in response to aluminum was revealed in maize and rice; these TFs were involved in phytohormone signaling and growth regulation [62][63][64]. In our study, expression of transcript TR35721|c0 g1, which encodes Agamous-like MADS-box protein AGL62, was induced by Al only in resistant flax genotypes (the binary logarithm of fold change, log , was equal to 2.3 and 3.0 in Hermes and TMP1919, respectively). The same was true for transcript TR25300|c0 g1 encoding NAC domain-containing protein 100: significant upregulation was revealed under aluminum stress in resistant cultivars (log was equal to 1.0 and 1.1 in Hermes and TMP1919, respectively), while some decrease of the expression-in sensitive ones. Thus, overexpression of genes encoding NAC domain-containing protein 100 and Agamous-like MADS-box protein AGL62 only in resistant flax genotypes under Al 3+ exposure indicates that these TFs could contribute to flax tolerance to aluminum. We also analyzed mRNA level alterations of other TFs, whose role in response to Al is known. The role of WRKY TFs in Al tolerance in plants was shown to be ambiguous: in Arabidopsis, WRKY46 was identified as a negative regulator of ALMT1 and mutations in WRKY46 lead to increased malate secretion and higher Al resistance [65]; on the contrary, WRKY22 promoted aluminum tolerance in rice through activation of FRDL4 (Ferric reductase defective 4) and increased citrate secretion [66]. In flax, we revealed slight upregulation of a number of WRKY-encoding genes but did not observe differences between resistant and sensitive cultivars. For genes encoding sensitive to proton rhizotoxicity 1 (STOP1) and Abscisic acid stress ripening (ASR), which are also known TFs associated with aluminum response [17,20,[67][68][69], no significant expression alterations were revealed by us in flax under aluminum stress. Thus, it can be suggested that Agamous-like MADS-box protein AGL62 and NAC domain-containing protein 100 are the leading TFs involved in Al tolerance in flax via regulation of plant growth and development.
Besides TFs, we observed more pronounced upregulation of genes encoding polygalacturonase, cellulose synthase, pectinesterase, beta-glucosidase, and GDSL esterase in resistant to aluminum flax cultivars compared to sensitive ones.
These enzymes are involved in cell wall metabolism, which plays important role in plant response to heavy metals and other abiotic stresses [70]. Modification of cell wall and changes in binding properties of the apoplast are known to contribute to Al tolerance in plants [70,71]. In flax, upregulation of genes encoding cell wall-related proteins was observed in resistant genotypes; therefore, cell wall modification could be one of the mechanisms of flax tolerance to Al ions.
In the top of the list of transcripts that were differentially expressed between resistant and sensitive genotypes (S8 Table), we also observed peroxidase and ABC transporterencoding genes, whose role in plant tolerance to aluminum stress is known [20]. Peroxidases are involved in detoxification of ROS and their overexpression in some plant species is associated with Al tolerance [20]. Moreover, overexpression of Arabidopsis peroxidase in tobacco plants improved their tolerance to aluminum [72]. At the same time, peroxidase expression was decreased under aluminum treatment in Camellia sinensis [73]. In wheat, peroxidase expression was induced by Al stress; however, their activity was lower in resistant genotypes compared to sensitive ones [74,75], while in chickpea, peroxidase activity was almost similar in tolerant and sensitive genotypes [76]. In our study on flax, for TR33816|c0 g1 transcript encoding peroxidase 5, a significant expression decrease was revealed in resistant genotypes (log was equal to -1.7 for both Hermes and TMP1919), while slight decrease or retention was identified in sensitive ones (log = −0.2 for Lira and log = −0.5 for Orshanskiy). The role of ABC transporters in plant response to aluminum stress is well characterized [8,[77][78][79][80][81]. In flax, ABC transporter-encoding transcript TR4576|c0 g2 had diverse expression changes between resistant and sensitive to aluminum genotypes: strong expression decrease was revealed in resistant cultivars (log was equal to -1.7 and -1.2 for Hermes and TMP1919, respectively) and slight expression decrease was observed in sensitive ones (log = −0.3 for Lira and log = −0.4 for Orshanskiy). For ABC transporter, we did not find associations of flax tolerance to Al with high expression levels of transcripts both under stress and control conditions. Considering the fact that expression of aluminum tolerance genes is usually higher in resistant genotypes and often increased by Al treatment [8], we suggested that peroxidase-and ABC transporterencoding genes are not the aluminum tolerance gene in flax.
The scoring for identification of differentially expressed genes in resistant genotypes compared to sensitive ones under Al 3+ exposure was shown to be the promising tool for revelation flax genes that are involved in tolerance to aluminum. Obtained results allowed us to suggest that some genes (including STOP1-, ASR-, peroxidase-, and ABC transporterencoding ones), whose role in aluminum response was revealed for several plant species, do not play the key role in flax tolerance to the stress. At the same time, MADS-box and NAC TFs, which regulate plant growth and development, and the enzymes that are involved in cell wall modifications are likely important for flax tolerance to aluminum. (1) (1)

Control conditions
Aluminum exposure

GO Analysis and Gene Expression Profiles.
For deeper understanding of the mechanisms of flax response to aluminum, we performed GO enrichment analysis to identify the processes in which genes with significant expression alterations under Al 3+ exposure are involved. Overrepresented GO terms were assessed for top lists (tops) of upand downregulated genes in flax plants under aluminum exposure. Tops of 50, 100, 300, and more than 300 differentially expressed genes (maximum number of genes was limited by the statistical significance of expression alterations corresponding to p-value < 0.05 in S1, S2, and S3 Tables) were used in the analysis. For the pool of all studied cultivars, overrepresented GO terms for the tops of 50 and 100 upregulated genes were related to circadian rhythm, multicellular organismal movement, and potassium: sodium symporter activity, while for the top of 300 upregulated genes, transferase activity had also appeared (S9 Table). For all the tops of downregulated genes for the pool of all studied genotypes, GO terms related to channel activity and transport were the most abundant. Heatmap for "water transmembrane transporter activity" GO term is presented in Figure 2 as an example. In the pool of sensitive to aluminum flax genotypes, we identified less number of significantly overrepresented GO terms compared to the pool of resistant ones for the tops of upregulated genes (S10 and S11 Tables). These GO terms were associated with apoplast, rhythmic processes, peptidase and oxidoreductase activities. In resistant cultivars, overrepresented GO terms included oxidoreductase activity and rhythmic processes too, but unlike sensitive genotypes, GO terms related to glutathione metabolism, glucosyltransferase activity, and response to stimulus were also overrepresented.
We also performed GO analysis for individual cultivars to evaluate the similarity in response to aluminum of genotypes of the same tolerance group (S12, S13, S14, and S15 Tables). For the top of 50 upregulated genes, glutathione transferase activity term was overrepresented in both resistant cultivars . Each heatmap row corresponds to one gene (total 21 000 genes). Color scale represents the binary logarithm of expression level fold change (aluminum exposure/control conditions) from -2 (i.e., 4-fold downregulation, blue) to +2 (4-fold upregulation, red).
(Hermes and TMP1919). However, for TMP1919, unlike Hermes, oxidoreductase activity term was also overrepresented and, in this, it was similar to sensitive cultivar Orshanskiy. When the top of 300 upregulated genes was used for analysis, much more GO terms were overrepresented. Glutathione transferase activity was in the top of the list not only in resistant genotypes but also in sensitive cultivar Orshanskiy, while oxidoreductase activity term was overrepresented in all studied genotypes. Gene expression profiles of particular GO terms had no relation to the degree of tolerance to aluminum. For example, expression alterations of genes related to oxidoreductase activity were more similar for pairs of cultivars TMP1919/Lira and Hermes/Orshanskiy (Figure 3).
GO enrichment analysis allowed us to reveal biological processes that are the most affected by aluminum stress in flax. Summarizing the results of this analysis, it could be concluded that Al 3+ exposure had negative effects on transmembrane transport in flax. It is not surprising since aluminum disturbs ion transport in plants [19]. It is known that oxidase overexpression confers aluminum tolerance of plants [82,83] and genes related to oxidoreductase activity GO term were upregulated in response to aluminum in flax. For glutathione transferase activity GO term, alterations in gene expression were more pronounced in resistant cultivars. Their role in Al response was discussed in our previous study [42].
To reveal coexpression gene networks in flax under aluminum stress, weighted gene coexpression network analysis (WGCNA) was performed for 8986 genes with CPM more than 20 in 6 or more samples. We identified 23 clusters (modules) of coexpressed genes based on similarity in gene expression (S16 Figure). Modules included from 46 to 1831 genes: module 1 was formed by 1831 genes enriched for RNA metabolic process and proton antiporter activity; module 2-organelle membrane and oxidationreduction process; 3-ribosome and macromolecule biosynthetic process; 6-transporter activity and rhythmic process; 7-regulation of cellular process and macromolecule biosynthesis; 8-phosphorus metabolic process, ion binding, and oxidation-reduction process; 9-nucleolus, ribosome, and RNA metabolic process; 13-mitochondrial part; 14-ion transport; 16-lipid metabolic processes. Genes from the same module are probably coregulated, and the processes, in which they are involved, may be related to each other. For other modules, we did not reveal significant enrichment in genes of particular GO terms.

Conclusions
We analyzed transcriptomes of flax cultivars with diverse tolerance to aluminum under control conditions and Al stress and identified genes that were significantly up-or downregulated. A number of genes whose expression alterations differed between resistant and sensitive cultivars were revealed, including those encoding MADS-box and NAC TFs and cell wall biogenesis related enzymes, suggesting them to be involved in aluminum tolerance in flax. These results indicate that, for flax, the internal tolerance to aluminum via reduction of Al ions toxicity is inherited. GO enrichment analysis led to the identification of processes that are affected by aluminum in both resistant and sensitive flax genotypes. The genotype effect on the response to Al 3+ exposure was strong enough; although resistant and sensitive cultivars had diverse gene expression changes under the stress, significant impact of a particular genotype on stress response was also observed. Understanding the mechanisms of flax response to Al and identification of tolerance genes is the basis for the development of improved cultivars which will retain high productivity on acid soils that is more preferable strategy than liming the soils.

Data Availability
Sequencing data used to support the findings of this study have been deposited in Sequence Read Archive-SRP089959.

Conflicts of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as potential conflicts of interest.
Supplementary Materials S1 Table. Gene expression alterations in all analyzed flax genotypes under aluminum stress. S2 Table. Gene expression alterations in resistant flax genotypes under aluminum stress. S3 Table. Gene expression alterations in sensitive flax genotypes under aluminum stress. S4 Table. Gene expression alterations in resistant flax cultivar Hermes under aluminum stress. S5 Table. Gene expression alterations in resistant flax cultivar TMP1919 under aluminum stress. S6 Table. Gene expression alterations in sensitive flax cultivar Lira under aluminum stress. S7 Table. Gene expression alterations in sensitive flax cultivar Orshanskiy under aluminum stress. S8 Table. Diverse expression alterations between groups of resistant and sensitive to aluminum flax genotypes. S9 Table. GO enrichment analysis for tops of up-and downregulated genes in the pool of all studied flax cultivars after aluminum exposure. S10 Table. GO enrichment analysis for tops of up-and downregulated genes in the pool of resistant flax cultivars after aluminum exposure. S11 Table. GO enrichment analysis for tops of up-and downregulated genes in the pool of sensitive flax cultivars after aluminum exposure. S12 Table. GO enrichment analysis for tops of up-and downregulated genes in resistant flax cultivar Hermes after aluminum exposure. S13 Table. GO enrichment analysis for tops of up-and downregulated genes in resistant flax cultivar TMP1919 after aluminum exposure. S14 Table. GO enrichment analysis for tops of up-and downregulated genes in sensitive flax cultivar Lira after aluminum exposure. S15 Table. GO enrichment analysis for tops of up-and downregulated genes in sensitive flax cultivar Orshanskiy after aluminum exposure. S16 Figure.