Transcriptomic Analysis for Indica and Japonica Rice Varieties under Aluminum Toxicity

Aluminum (Al) at high concentrations inhibits root growth, damage root systems, and causes significant reductions in rice yields. Indica and Japonica rice have been cultivated in distinctly different ecological environments with different soil acidity levels; thus, they might have different mechanisms of Al-tolerance. In the present study, transcriptomic analysis in the root apex for Al-tolerance in the seedling stage was carried out within Al-tolerant and -sensitive varieties belonging to different subpopulations (i.e., Indica, Japonica, and mixed). We found that there were significant differences between the gene expression patterns of Indica Al-tolerant and Japonica Al-tolerant varieties, while the gene expression patterns of the Al-tolerant varieties in the mixed subgroup, which was inclined to Japonica, were similar to the Al-tolerant varieties in Japonica. Moreover, after further GO (gene ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) analyses of the transcriptomic data, we found that eight pathways, i.e., “Terpenoid backbone biosynthesis”, “Ribosome”, “Amino sugar and nucleotide sugar metabolism”, “Plant hormone signal transduction”, “TCA cycle”, “Synthesis and degradation of ketone bodies”, and “Butanoate metabolism” were found uniquely for Indica Al-tolerant varieties, while only one pathway (i.e., “Sulfur metabolism”) was found uniquely for Japonica Al-tolerant varieties. For Al-sensitive varieties, one identical pathway was found, both in Indica and Japonica. Three pathways were found uniquely in “Starch and sucrose metabolism”, “Metabolic pathway”, and “Amino sugar and nucleotide sugar metabolism”.


Introduction
Over 50% of the world's arable land is acidic, and about 13% of global rice is produced on acidic soils [1]. The acidity of some soils is due to aluminum (Al), which is the most highly abundant metal in the earth's crust. Acidic soil containing solubilized trivalent Al (Al 3+ ) will have pH values lower than 5.0. Not only is the growth of rice roots inhibited, but rice root systems can also be damaged by high concentrations of Al in the soil, which can both lead to significant reductions in rice yields [2,3]. Uncovering the genetic mechanisms of Al-tolerance in rice is the premise for rice Al-tolerance improvement.

The Identification of Varieties with Extreme Al-Tolerances and Susceptibilities, According to the Population Structure in Ting's Core Collection
The landraces with extreme Al-tolerance and susceptibility in each subgroup in Ting's core collection were selected by measuring the relative root elongation (RRE) in the presence of Al toxicity (Figure 1), i.e., Ba shi zi (Al-tolerant) and Chang ning wu qu nan tou zhan (Al-sensitive) in the subgroup of Indica, Ai you (Al-tolerant) and Kai xuan (Al-sensitive) in the subgroup of Japonica, and Bnlastog (Al-tolerant) and Hei ke da nuo (Al-sensitive) in the mixed subgroup.

Transcriptome Sequencing and Sequence Quality
There were 430,562,708 raw reads gained in total. The following quality control steps were performed: adapters were cut from the reads; reads having N base read frequencies of greater than 10% were filtered out; 3 -end base pairs with Q < 30 were cut; and low-quality reads were filtered out. The remaining 429,380,993 reads were used for downstream analyses (Table 1).
We performed a correlation analysis for the three biological replications, and the correlation coefficients ranged from 0.792 to 0.945 (p < 0.01, Figure 2), which indicates a high correlation among the three biological replications. In this study, we used fragments per kilobase of the exon model per million mapped reads (FPKM), to demonstrate the expression levels of the genes. We defined a gene to be significantly differentially expressed when a gene's log 2 (FPKM after Al treatment/FPKM without Al treatment) was higher than +1 and lower than −1. The positive value of log 2 (FPKM_after_Al treatment/FPKM_before_Al) showed that the expression of the gene was up-regulated, while the negative value showed that it was down-regulated.

Transcriptome Sequencing and Sequence Quality
There were 430,562,708 raw reads gained in total. The following quality control steps were performed: adapters were cut from the reads; reads having N base read frequencies of greater than 10% were filtered out; 3′-end base pairs with Q < 30 were cut; and low-quality reads were filtered out. The remaining 429,380,993 reads were used for downstream analyses (Table 1). We performed a correlation analysis for the three biological replications, and the correlation coefficients ranged from 0.792 to 0.945 (p < 0.01, Figure 2), which indicates a high correlation among the three biological replications. In this study, we used fragments per kilobase of the exon model per million mapped reads (FPKM), to demonstrate the expression levels of the genes. We defined a gene to be significantly differentially expressed when a gene's log2(FPKM after Al treatment/FPKM without Al treatment) was higher than +1 and lower than −1. The positive value of log2(FPKM_after_Al treatment/FPKM_before_Al) showed that the expression of the gene was upregulated, while the negative value showed that it was down-regulated.

Comparison of Significantly Expressed Genes among the Three Al-Tolerant Varieties
There were 3301 genes that were significantly up-regulated, both in Ba shi zi and Ai you, 4837 genes significantly up-regulated both in Ba shi zi and Bnlastog, 2695 genes significantly up-regulated both in Ai you and Bnlastog, and 2484 genes significantly up-regulated in all three varieties. There were 1294 genes that were uniquely up-regulated in Ba shi zi, 647 genes uniquely up-regulated in Ai

Comparison of Significantly Expressed Genes among the Three Al-Tolerant Varieties
There were 3301 genes that were significantly up-regulated, both in Ba shi zi and Ai you, 4837 genes significantly up-regulated both in Ba shi zi and Bnlastog, 2695 genes significantly up-regulated both in Ai you and Bnlastog, and 2484 genes significantly up-regulated in all three varieties. There were 1294 genes that were uniquely up-regulated in Ba shi zi, 647 genes uniquely up-regulated in Ai you, and 1550 genes uniquely up-regulated in Bnlastog. There were 854 genes significantly down-regulated both in Ba shi zi and Ai you, 3309 genes significantly down-regulated both in Ba shi zi and Bnlastog, 768 genes significantly down-regulated both in Ai you and Bnlastog, and 709 genes significantly down-regulated in all three varieties. There were 1456 genes uniquely down-regulated in Ba shi zi, 257 genes uniquely up-regulated in Ai you, and 973 genes uniquely down-regulated in Bnlastog ( Figure 3A).

Gene Expression Patterns of Al-Tolerant and Al-Sensitive Varieties
In order to further analyze the differences in Al-tolerance between Indica and Japonica, gene expression patterns of Al-tolerant and Al-sensitive varieties in three subgroups under Al toxicity were identified. Ba shi zi of Indica Al-tolerant and Ai you of Japonica Al-tolerant had significantly different expression patterns before and after applying the Al stress, and Bnlastog of mixed Altolerant had a similar expression pattern to Ba shi zi ( Figure 4A). According to Cheng's index, Bnlastog is more closely related to Indica (Table 2). In addition, Chang ning wu qu nan tou zhan of Indica Al-sensitive and Kai xuan of Japonica Al-sensitive had significantly different expression patterns before and after applying Al stress, and Hei ke da nuo of mixed Al-sensitive had a similar expression pattern to Kai xuan ( Figure 4B). According to Cheng's index, Hei ke da nuo is more closely related to Japonica (Table 2).

Comparison of Significantly Expressed Genes among the Three Al-Sensitive Varieties
There were 3012 genes that were significantly up-regulated, both in Chang ning wu qu nan tou zhan and Kai xuan, 710 genes significantly up-regulated both in Chang ning wu qu nan tou zhan and Hei ke da nuo, 902 genes significantly up-regulated both in Kai xuan and Hei ke da nuo, and 682 genes significantly up-regulated in all three varieties. There were 1659 genes uniquely up-regulated in Chang ning wu qu nan tou zhan, 1338 genes uniquely up-regulated in Kai xuan, and 133 genes uniquely up-regulated in Hei ke da nuo. There were 1410 genes significantly down-regulated both in Chang ning wu qu nan tou zhan and Kai xuan, 258 genes significantly down-regulated both in Chang ning wu qu nan tou zhan and Hei ke da nuo, 374 genes significantly down-regulated both in Kai xuan and Hei ke da nuo, and 709 genes significantly down-regulated in all three varieties. There were 1713 genes uniquely down-regulated in Chang ning wu qu nan tou zhan, 396 genes uniquely down-regulated in Kai xuan, and 150 genes uniquely down-regulated in Hei ke da nuo ( Figure 3B).

Gene Expression Patterns of Al-Tolerant and Al-Sensitive Varieties
In order to further analyze the differences in Al-tolerance between Indica and Japonica, gene expression patterns of Al-tolerant and Al-sensitive varieties in three subgroups under Al toxicity were identified. Ba shi zi of Indica Al-tolerant and Ai you of Japonica Al-tolerant had significantly different expression patterns before and after applying the Al stress, and Bnlastog of mixed Al-tolerant had a similar expression pattern to Ba shi zi ( Figure 4A). According to Cheng's index, Bnlastog is more closely related to Indica (Table 2). In addition, Chang ning wu qu nan tou zhan of Indica Al-sensitive and Kai xuan of Japonica Al-sensitive had significantly different expression patterns before and after applying Al stress, and Hei ke da nuo of mixed Al-sensitive had a similar expression pattern to Kai xuan ( Figure 4B). According to Cheng's index, Hei ke da nuo is more closely related to Japonica (Table 2). expression patterns of Al-tolerant and Al-sensitive varieties in three subgroups under Al toxicity were identified. Ba shi zi of Indica Al-tolerant and Ai you of Japonica Al-tolerant had significantly different expression patterns before and after applying the Al stress, and Bnlastog of mixed Altolerant had a similar expression pattern to Ba shi zi ( Figure 4A). According to Cheng's index, Bnlastog is more closely related to Indica (Table 2). In addition, Chang ning wu qu nan tou zhan of Indica Al-sensitive and Kai xuan of Japonica Al-sensitive had significantly different expression patterns before and after applying Al stress, and Hei ke da nuo of mixed Al-sensitive had a similar expression pattern to Kai xuan ( Figure 4B). According to Cheng's index, Hei ke da nuo is more closely related to Japonica (Table 2).  Note: The Cheng's index is between 1-7, 8-13, 14-17, and 18-24 for typical Indica, Indica-clined, Japonica-clined, and typical Japonica rice, respectively.  Note: The Cheng's index is between 1-7, 8-13, 14-17, and 18-24 for typical Indica, Indica-clined, Japonica-clined, and typical Japonica rice, respectively.

Classification of Differentially Expressed Genes of Al-Tolerant Varieties in Indica and Japonica
As shown with the gene ontology (GO) analyses of the transcriptome data, Ba shi zi, in comparison with the Al-tolerant variety Ai you, had more differentially expressed genes categorized as a "cellular component", but fewer differentially expressed genes categorized as "molecular function" and "biological process" ( Figure 5A). In addition, Chang ning wu qu nan tou zhan, in comparison with Kai xuan, had more differentially expressed genes categorized as "cellular component", but fewer differentially expressed genes categorized as "molecular function" and "biological process" ( Figure 5B). Besides, pathway enrichment analysis of differentially expressed genes was performed using g:Profiler. We found that genes were uniquely enriched for Bas hi zi into "structural molecule activity (GO:0005198)", "transferase activity, transferring hexosyl groups (GO:0016758)", "structural constituent of ribosome (GO:0003735)", "transferase activity, transferring hexosyl groups (GO:0016757)" and "catalytic activity (GO:0003824)" in the category for "Molecular function", as well as "Carbohydrate metabolic process (GO:0005975)" in the category for "Biological process". Genes were uniquely enriched for Ai you and Bas hi zi into the "cell periphery (GO:0071944)", and "ribosome (GO:0005840)" in the category of "cellular component", respectively (Table S1).
(GO:0016758)", "structural constituent of ribosome (GO:0003735)", "transferase activity, transferring hexosyl groups (GO:0016757)" and "catalytic activity (GO:0003824)" in the category for "Molecular function", as well as "Carbohydrate metabolic process (GO:0005975)" in the category for "Biological process". Genes were uniquely enriched for Ai you and Bas hi zi into the "cell periphery (GO:0071944)", and "ribosome (GO:0005840)" in the category of "cellular component", respectively (Table S1). Moreover, we conducted pathway-based analysis for differentially expressed genes by the KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway database. Six pathways were detected, both in Indica and Japonica Al-tolerant varieties. Eight pathways, i.e., "Terpenoid backbone biosynthesis", "Ribosome", "Amino sugar and nucleotide sugar metabolism", "Plant hormone signal transduction", "TCA cycle", "Synthesis and degradation of ketone bodies", and "Butanoate metabolism" were found uniquely for the Indica Al-tolerant variety Ba shi zi, while only one pathway (i.e., "Sulfur metabolism") was found uniquely for the Japonica Al-tolerant variety Ai you ( Figure  6A,B). In the Al-sensitive sub-species, one identical pathway was found, both in Indica and Japonica. Three pathways were found uniquely in Chang ning wu qu nan tou zhan, i.e., "Starch and sucrose metabolism", "Metabolic pathway", and "Amino sugar and nucleotide sugar metabolism". Four Alsensitive pathways were detected in Indica, and one was detected in Japonica. Hence, there was no figure to show the enrichment of the KEGG pathway for these two Al-sensitive varieties in the present study. The detailed information of the above KEGG pathways is summarized in Tables S2-S4. Moreover, we conducted pathway-based analysis for differentially expressed genes by the KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway database. Six pathways were detected, both in Indica and Japonica Al-tolerant varieties. Eight pathways, i.e., "Terpenoid backbone biosynthesis", "Ribosome", "Amino sugar and nucleotide sugar metabolism", "Plant hormone signal transduction", "TCA cycle", "Synthesis and degradation of ketone bodies", and "Butanoate metabolism" were found uniquely for the Indica Al-tolerant variety Ba shi zi, while only one pathway (i.e., "Sulfur metabolism") was found uniquely for the Japonica Al-tolerant variety Ai you ( Figure 6A,B). In the Al-sensitive sub-species, one identical pathway was found, both in Indica and Japonica. Three pathways were found uniquely in Chang ning wu qu nan tou zhan, i.e., "Starch and sucrose metabolism", "Metabolic pathway", and "Amino sugar and nucleotide sugar metabolism". Four Al-sensitive pathways were detected in Indica, and one was detected in Japonica. Hence, there was no figure to show the enrichment of the KEGG pathway for these two Al-sensitive varieties in the present study. The detailed information of the above KEGG pathways is summarized in Tables S2-S4. Besides enrichment analysis based on the KEGG pathway, further analysis was performed with the transcriptomic data, and 45 genes were found that had differential expressions between the Indica and Japonica subgroups, which merit further study. These included nine genes relating to abiotic stress, four relating to stress, four relating to the binding of heavy metals, 20 relating to cell walls, and three relating to gibberellin (Table 3). Besides enrichment analysis based on the KEGG pathway, further analysis was performed with the transcriptomic data, and 45 genes were found that had differential expressions between the Indica and Japonica subgroups, which merit further study. These included nine genes relating to abiotic stress, four relating to stress, four relating to the binding of heavy metals, 20 relating to cell walls, and three relating to gibberellin (Table 3). Silicon transporter involved in the distribution of silicon in shoots. Is responsible for the transport of silicon from the xylem to the leaf tissues. Silicon is beneficial to plant growth and helps plants to overcome abiotic and biotic stresses by preventing lodging (falling over), and increasing plant resistance to pests and diseases, as well as other stresses (PubMed:18515498). In the nodes, is involved with LSI2 and LSI3 in silicon intervascular transfer, which is required for the preferential distribution of silicon, such as the hyperaccumulation of silicon in the husk

LOC_Os02g44230
Abiotic stress 3.78 no Removes the phosphate from trehalose 6-phosphate to produce free trehalose. Trehalose accumulation in the plant improves abiotic stress tolerance  Thought to be a Golgi-localized beta-glycan synthase that polymerize the backbones of noncellulosic polysaccharides (hemicelluloses) of the plant cell wall. Required for the synthesis of a cell wall polysaccharide that is essential for root hair elongation, but not initiation. May be the functional ortholog of Arabidopsis CSLD3/KOJAK

LOC_Os01g60770
Cell wall −0.83 1.10 May cause loosening and extension of plant cell walls by disrupting non-covalent bonding between cellulose microfibrils and matrix glucans. No enzymatic activity has been found. May be required for rapid internodal elongation in deepwater rice during submergence (By similarity)

LOC_Os05g19570
Cell wall −1.56 1.30 May cause loosening and extension of plant cell walls by disrupting non-covalent bonding between cellulose microfibrils and matrix glucans. No enzymatic activity has been found. May be required for rapid internodal elongation in deepwater rice during submergence (By similarity)

LOC_Os05g19600
Cell wall −4.13 1.79 May cause loosening and extension of plant cell walls by disrupting non-covalent bonding between cellulose microfibrils and matrix glucans. No enzymatic activity has been found. May be required for rapid internodal elongation in deepwater rice during submergence (By similarity)

LOC_Os03g05110
Cell wall −0.60 0.92 Involved in the attachment of the Gal residue on the third xylosyl unit within the XXXG core structure of xyloglucan, the principal glycan that interlaces the cellulose microfibrils in plant cell wall. Interacts with actin, and is required for proper endomembrane organization, and for the cell elongation (By similarity)

LOC_Os05g43530
Cell wall −2.06 no Probable beta-1,4-glucan synthase, involved in the synthesis of the xyloglucan backbone instead of cellulose. Seems to work simultaneously with xyloglucan 6-xylosyltransferase. Xyloglucan is a type of noncellulosic polysaccharide g of the plant cell wall, and it consists of a glucan backbone substituted by xylose, galactose, and fucose (By similarity) May cause loosening and extension of plant cell walls by disrupting non-covalent bonding between cellulose microfibrils and matrix glucans. No enzymatic activity has been found. May be required for rapid internodal elongation in deepwater rice during submergence (By similarity)

LOC_Os04g46650
Cell wall −3.05 no May cause loosening and extension of plant cell walls by disrupting non-covalent bonding between cellulose microfibrils and matrix glucans. No enzymatic activity has been found. May be required for rapid internodal elongation in deepwater rice during submergence (By similarity)

LOC_Os03g44290
Cell wall −3.20 no May cause loosening and extension of plant cell walls by disrupting non-covalent bonding between cellulose microfibrils and matrix glucans. No enzymatic activity has been found. May be required for rapid internodal elongation in deepwater rice during submergence (By similarity)

LOC_Os05g15690
Cell wall May cause loosening and extension of plant cell walls by disrupting non-covalent bonding between cellulose microfibrils and matrix glucans. No enzymatic activity has been found. May be required for rapid internodal elongation in deepwater rice during submergence (By similarity) Catalyzes the 3-beta-hydroxylation of the inactive gibberellin precursors, leading to the formation of bioactive gibberellins. In vitro, converts the precursors GA20, GA5, GA44, and GA9 to the corresponding 3-beta-hydroxylated active products GA1, GA3, GA38, and GA4, respectively. Involved in the production of bioactive GA for vegetative growth and development (PubMed:11438692). Controls the elongation of the vegetative shoot and plant height by the regulation of active gibberellin levels

LOC_Os05g06670
Gibberellin 2.53 no Catalyzes the 2-beta-hydroxylation of several biologically active gibberellins, leading to the homeostatic regulation of their endogenous levels. Catabolism of gibberellins (GAs) plays a central role in plant development. Controls the levels of bioactive GAs in the shoot apical meristem, which regulate the vegetative-to-reproductive phase transition. In vitro, converts GA1, GA4, GA9, GA20, and GA44 to the corresponding 2-beta-hydroxylated products GA8, GA34, GA51, GA29, and GA98, respectively Catalyzes three successive oxidations of the 4-methyl group of ent-kaurene, giving kaurenoic acid, a key step in gibberellins (GAs) biosynthesis. GAs, which are involved in many processes, including stem elongation, play a central role in plant development

qRT-PCR Verification
We randomly selected eight genes with significantly differential expressions after applying Al stress in Al-tolerant and Al-sensitive varieties. The expression levels of the above eight genes were verified by qRT-PCR (Quantitative reverse transcription-PCR) analysis, and they were significantly different in Nipponbare and Kasalath (Figure 7), which showed that the transcriptome sequencing results were reliable.

qRT-PCR Verification
We randomly selected eight genes with significantly differential expressions after applying Al stress in Al-tolerant and Al-sensitive varieties. The expression levels of the above eight genes were verified by qRT-PCR (Quantitative reverse transcription-PCR) analysis, and they were significantly different in Nipponbare and Kasalath (Figure 7), which showed that the transcriptome sequencing results were reliable.

Study of Al-Tolerance within Varieties in Different Subgroups, Based on Population Structure
Al-tolerant and -sensitive accessions were discovered, both in Indica and Japonica rice [7,10]. Moreover, Al-tolerant mechanisms of rice were uncovered within Indica and Japonica rice in previous studies. For instance, one Al-tolerant Japonica variety named Azucena, and one Al-sensitive Indica variety named IR1552 were used in mapping Al-related QTL in rice [15,19]. Other studies focused on Al-tolerant and Al-sensitive Indica varieties named Xiangnuo 1 and Xiangzhongxian 2, respectively [20][21][22]. Huang et al. and Yamaji et al. uncovered Al-tolerance mechanisms in a Japonica variety, originating from Koshihikari, using one Al-sensitive and one Al-tolerant mutant, respectively [4,5]. One Al-tolerant Japonica variety, named Nipponbare, was used to identify Al-tolerant mechanisms in rice [23]. Wagatsuma et al. investigated the Al-tolerance mechanism with five temperate Japonica varieties, including two Al-sensitive and three Al-tolerant varieties, as well as one Al-sensitive Indica variety [24]. However, the Al-tolerant genotypes used in the previous studies were Japonica or mutants originating from Japonica. Besides, we speculated that the mechanisms between Japonica and Indica might not be common, according to the distinct ecological environment, especially the acidity level in the soil for Indica and Japonica rice. Thus, Al-tolerant and Al-sensitive varieties in the different subgroups, according to the population structure of Ting's core collection [25], were used to dissect the mechanisms of Al-tolerance in rice, in the present study. Choosing Indica-or Japonica-specific genes relating to Al-tolerance, using marker-assisted selection (MAS), would be very significant for the improvement of rice Al-tolerance.

Different Gene Expression Patterns of Al-Tolerance between Indica and Japonica
After the transcriptomic analysis, different gene expression patterns of Al-tolerant and Alsensitive varieties in Indica and Japonica under Al toxicity were uncovered. Indica Al-tolerant and Japonica Al-tolerant strains had significantly different expression patterns before and after applying Al stress, while mixed Al-tolerant, which was inclined to Indica, had a similar expression pattern to Indica Al-tolerant. Moreover, Indica Al-sensitive and Japonica Al-sensitive strains had significantly different expression patterns before and after applying Al stress, and the mixed Al-sensitive, which Figure 7. Expression of eight candidate genes in the rice root tip. * and ** represent significant differences at the p < 0.05 and p < 0.01 levels (t-test), respectively.

Study of Al-Tolerance within Varieties in Different Subgroups, Based on Population Structure
Al-tolerant and -sensitive accessions were discovered, both in Indica and Japonica rice [7,10]. Moreover, Al-tolerant mechanisms of rice were uncovered within Indica and Japonica rice in previous studies. For instance, one Al-tolerant Japonica variety named Azucena, and one Al-sensitive Indica variety named IR1552 were used in mapping Al-related QTL in rice [15,19]. Other studies focused on Al-tolerant and Al-sensitive Indica varieties named Xiangnuo 1 and Xiangzhongxian 2, respectively [20][21][22]. Huang et al. and Yamaji et al. uncovered Al-tolerance mechanisms in a Japonica variety, originating from Koshihikari, using one Al-sensitive and one Al-tolerant mutant, respectively [4,5]. One Al-tolerant Japonica variety, named Nipponbare, was used to identify Al-tolerant mechanisms in rice [23]. Wagatsuma et al. investigated the Al-tolerance mechanism with five temperate Japonica varieties, including two Al-sensitive and three Al-tolerant varieties, as well as one Al-sensitive Indica variety [24]. However, the Al-tolerant genotypes used in the previous studies were Japonica or mutants originating from Japonica. Besides, we speculated that the mechanisms between Japonica and Indica might not be common, according to the distinct ecological environment, especially the acidity level in the soil for Indica and Japonica rice. Thus, Al-tolerant and Al-sensitive varieties in the different subgroups, according to the population structure of Ting's core collection [25], were used to dissect the mechanisms of Al-tolerance in rice, in the present study. Choosing Indicaor Japonica-specific genes relating to Al-tolerance, using marker-assisted selection (MAS), would be very significant for the improvement of rice Al-tolerance.

Different Gene Expression Patterns of Al-Tolerance between Indica and Japonica
After the transcriptomic analysis, different gene expression patterns of Al-tolerant and Al-sensitive varieties in Indica and Japonica under Al toxicity were uncovered. Indica Al-tolerant and Japonica Al-tolerant strains had significantly different expression patterns before and after applying Al stress, while mixed Al-tolerant, which was inclined to Indica, had a similar expression pattern to Indica Al-tolerant. Moreover, Indica Al-sensitive and Japonica Al-sensitive strains had significantly different expression patterns before and after applying Al stress, and the mixed Al-sensitive, which was inclined to Japonica, had a similar expression pattern as Al-sensitive Japonica ( Figure 4). As shown with a enrichment analysis of differentially expressed genes using website of g:Profiler (http://biit.cs.ut.ee/ gprofiler/), differentially expressed genes were uniquely enriched into different terms for Al-tolerant Indica and Japonica. Some of the terms were reported to be related to plant Al-tolerance. For instance, it was demonstrated that the over-expression of one Al-induced gene parB (glutathione S-transferase) could ameliorate Al toxicity in Arabidopsis [26]. Differentially expressed genes, which are enriched in the "Structural constituent of ribosome", were found to be related to Al-tolerance in Stylosanthes [27].

Different Pathways of Al-Tolerance between Indica and Japonica
Al-tolerance in plants is a complex process that is mediated by various factors. In the present study, 11 pathways, which might be the difference for Al-tolerance between Indica and Japonica rice, were detected. Some of the pathways were reported to relate to Al-tolerance in other species; it was reported that hormonal equilibrium could be regulated by nitric oxide to improve Al-tolerance in the root apices of rye and wheat [28]. Ninety-two genes related to the plant hormone signaling transduction, pyruvate metabolism, amino sugar and nucleotide sugar metabolism, ribosome, the TCA cycle, the synthesis and degradation of ketone bodies, and butanoate metabolism, were enriched uniquely in Al-tolerant Indica (Table S2) and most of them have been reported to play important roles in Al-tolerance. For instance, LOC_Os08g39830 and LOC_Os03g53150, which regulate ethylene-insensitivity, respectively, [29,30] were expressed differentially after Al treatment in our study, and ethylene, as well as auxin could mediate the Al-induced inhibition of root elongation in Arabidopsis [31]. LOC_Os10g33800 and LOC_Os10g33800, belonging to the KEGG pathway of pyruvate metabolism and the TCA cycle, respectively, regulate the synthesis of malate dehydrogenase, which was over-expressed to confer tolerance to Al in alfalfa (Medicago sativa) [32]. Three genes were enriched uniquely in the KEGG pathway of sulfur metabolism in Al-tolerant Japonica (Table S3) and this pathway has been reported to be related to Al-tolerance [33]. For Al-sensitive Indica, 26 genes related to the pathways of amino sugars and nucleotide sugar metabolism, metabolic pathways, and starch and sucrose metabolism were enriched uniquely. Regarding amino sugar and nucleotide sugar metabolism, it has been reported that sugar concentrations increased with Al toxic treatment in rice [34].

Plant Material
Six varieties were chosen for transcriptomic analysis, according to the population structure existing in Ting's core collection [25]: two varieties (Al-tolerant and -sensitive) in the Indica subpopulation, two varieties (Al-tolerant and -sensitive) in the Japonica subpopulation, and two varieties (Al-tolerant and -sensitive) in the mixed subpopulation.

Phenotyping for Al-Tolerance
All of the rice seeds were harvested from the farm of the China National Rice Research Institute, Hangzhou (30 • 3 N, 120 • 2 E), during the late season (May-October) in 2014. A 1% H 2 O 2 solution was used to surface-sterilize 50 uniform seed of each rice variety. Deionized water was then used to wash the seeds three times. Next, the fifty seeds were soaked in 30 • C temperature deionized water in a dark growth chamber for two days, in order for them to germinate. After germination, the seedlings were placed on a plastic net, which was put into a 1.5 L plastic container with 0.5 mmol·L −1 CaCl 2 (pH = 4.0) solution. We renew the CaCl 2 solution on a daily basis. Precipitation of Al and other elements in Yoshida's solution [35]. will complicate Al-tolerance screening. It has been proven in many studies that the CaCl 2 solution can effectively screen young seedlings whose seeds are still able to produce the essential mineral nutrients within Al-tolerant rice [6,9]. Thus, the CaCl 2 solution was used instead of Yoshida's solution in this study. After 48 hr, the most uniform 20 seedlings were selected and used for Al toxicity treatment. We expose the 20 selected seedlings to 0.5 mmol·L −1 CaCl 2 (pH = 4.0) with 100 µmol·L −1 AlCl 3 . We expose control seedlings to 0.5 mmol·L −1 CaCl 2 (pH = 4.0) only. We calculated RRE (RRE = root length under Al toxicity treatment/root length without Al toxicity treatment) after 24 hr, to evaluate the Al-tolerance of all varieties. Six replications were performed. During each replication, root lengths of 10 seedlings from each variety were measured using a ruler, and the measurements were recorded, both before the treatments and after the treatments (14 hr). To measure Al-tolerance, we calculate the RRE value for each accession of the Ting's core collection. The criterion for finding out Al-tolerant varieties was RRE ≥ 0.5 in the present study [36].

Transcriptomic Analysis
Total RNA of the root tip after 24 hr under 100 µmol·L −1 AlCl 3 was extracted with an miRNeasy Kit (QIAGEN, Germantown, MD, USA). The RNA samples were evaluated on agarose gels, quantified in a spectrophotometer, and stored at −80 • C. RNA sequencing was performed by Illumina NextSeq 500 (San Diego, CA, USA) for 2 × 125-bp paired-end sequencing. Samples under Al toxic treatment (three times), and without Al toxic treatment (one time) were collected for RNA sequencing.

Selecting Genes with Different Expression Patterns between Indica and Japonica
We defined genes with the following attributes as having different expression patterns among Indica, Japonica, and mixed: (1) Genes were significantly up-regulated in one subgroup, but significantly down-regulated in the other subgroup; (2) Genes were significantly up-regulated or down-regulated in one subgroup, with a log 2 (fold change) that was greater than or equal to 2, but it did not have significant differential expression in the other subgroup. Genes with different expression patterns among the three subgroups were annotated with gene functions downloaded from UniProt [37].

Gene Ontology Plots
We used WEGO 2.0 (http://wego.genomics.org.cn/) [38] to visualize the ontology of genes with different expression patterns between Indica and Japonica. GO plots for the tolerant varieties were plotted, using all genes with different expression levels between tolerant varieties in Indica and Japonica as candidates. In order for WEGO to generate a graph of genes with different expression levels between sensitive varieties in Indica and Japonica, genes were significantly up-regulated or down-regulated in one subgroup, with a log 2 (fold change) of less than 2, but those with no significant differential expression in the other subgroup were also used as candidates for WEGO.