Tolerance to mild salinity stress in japonica rice: A genome-wide association mapping study highlights calcium signaling and metabolism genes

Salinity tolerance is an important quality for European rice grown in river deltas. We evaluated the salinity tolerance of a panel of 235 temperate japonica rice accessions genotyped with 30,000 SNP markers. The panel was exposed to mild salt stress (50 mM NaCl; conductivity of 6 dS m-1) at the seedling stage. Eight different root and shoot growth parameters were measured for both the control and stressed treatments. The Na+ and K+ mass fractions of the stressed plants were measured using atomic absorption spectroscopy. The salt treatment affected plant growth, particularly the shoot parameters. The panel showed a wide range of Na+/K+ ratio and the temperate accessions were distributed over an increasing axis, from the most resistant to the most susceptible checks. We conducted a genome-wide association study on indices of stress response and ion mass fractions in the leaves using a classical mixed model controlling structure and kinship. A total of 27 QTLs validated by sub-sampling were identified. For indices of stress responses, we also used another model that focused on marker × treatment interactions and detected 50 QTLs, three of which were also identified using the classical method. We compared the positions of the significant QTLs to those of approximately 300 genes that play a role in rice salt tolerance. The positions of several QTLs were close to those of genes involved in calcium signaling and metabolism, while other QTLs were close to those of kinases. These results reveal the salinity tolerance of accessions with a temperate japonica background. Although the detected QTLs must be confirmed by other approaches, the number of associations linked to candidate genes involved in calcium-mediated ion homeostasis highlights pathways to explore in priority to understand the salinity tolerance of temperate rice.

Salinity tolerance is an important quality for European rice grown in river deltas. We evaluated the salinity tolerance of a panel of 235 temperate japonica rice accessions genotyped with 30,000 SNP markers. The panel was exposed to mild salt stress (50 mM NaCl; conductivity of 6 dS m -1 ) at the seedling stage. Eight different root and shoot growth parameters were measured for both the control and stressed treatments. The Na + and K + mass fractions of the stressed plants were measured using atomic absorption spectroscopy. The salt treatment affected plant growth, particularly the shoot parameters. The panel showed a wide range of Na + /K + ratio and the temperate accessions were distributed over an increasing axis, from the most resistant to the most susceptible checks. We conducted a genome-wide association study on indices of stress response and ion mass fractions in the leaves using a classical mixed model controlling structure and kinship. A total of 27 QTLs validated by subsampling were identified. For indices of stress responses, we also used another model that focused on marker × treatment interactions and detected 50 QTLs, three of which were also identified using the classical method. We compared the positions of the significant QTLs to those of approximately 300 genes that play a role in rice salt tolerance. The positions of several QTLs were close to those of genes involved in calcium signaling and metabolism, while other QTLs were close to those of kinases. These results reveal the salinity tolerance of accessions with a temperate japonica background. Although the detected QTLs must be confirmed by other approaches, the number of associations linked to candidate genes involved in calcium-mediated ion homeostasis highlights pathways to explore in priority to understand the salinity tolerance of temperate rice. PLOS  Introduction these mechanisms is vigorous growth to dilute salt concentrations in the tissues [22]. This tolerance mechanism is possibly present in Nona Bokra, which develops a huge biomass. Another mechanism to avoid rapid Na + accumulation is sodium exclusion from the transpiration stream, which occurs by sequestrating Na + into the older leaves [23]. Another mitigation mechanism is ion compartmentation into cell vacuoles, which is mediated by Na + /H + antiporters, such as those from the vacuolar Na+/H+ antiporter (OsNHX) family [24]. A final mechanism is scavenging of phytotoxic reactive oxygen species (ROS) produced under stress by the activation of antioxidant enzymes, such as ascorbate peroxidase or peroxide dismutase, which limits membrane injuries and protects chloroplast function [13]. Although some authors have indicated that none of these mechanisms are preferentially used in rice [18], other authors have suggested that tissue tolerance mechanisms play secondary roles [14]. When the traits that may be affected by the salinity stress are assessed irrespective of the involved mechanisms (i.e., in the whole-plant perspective), vegetative growth is reduced. Shoots are more affected than roots [13], which in turn increases the root-to-shoot ratio. Cell expansion in young leaves is affected, decreasing the specific leaf area [25]. The thickening of leaf cell walls is sometimes considered adaptive because this process provides a greater volume in which salt can be sequestered [26]. Specific leaf area linearly increases with increasing relative growth rate [25]. In addition, at a later stage, salt stress induces delayed panicle initiation and flowering. Yield components, notably spikelet fertility, are affected by Na + accumulation in flag leaves and developing panicles [26,27].
The genetic control of salt tolerance is only partially understood. Because of the diversity of the involved mechanisms, the genetic architecture of salt tolerance is complex. Numerous quantitative trait loci (QTLs) have been identified using bi-parental mapping populations (reviewed by Negrao et al. [13] for studies conducted prior to 2011, [28][29]), and, more recently, association mapping panels [12,[30][31][32][33]. Lists of genes involved in salinity tolerance based on the study of mutants or transgenic plants have been established [13,34]. The gene OsHKT1;5, also called SKC1, which encodes a Na + -selective transporter, has been implicated in the tolerance of Nona Bokra [35]. A major QTL from Pokkali, named Saltol, which is involved in Na/K homeostasis, has been fine mapped in the same area of chromosome 1, but genes from this segment other than OsHKT1;5 may be implicated as well [36][37]. Pokkali carries the same haplotype as Nona Bokra at HKT1;5 [38]. The tolerance allele of HKT1;5 may originate from the aromatic groups [14]. Pokkali and Nona Bokra have been largely used as donors in breeding programs for the marker-assisted selection of salinity tolerance and varieties as tolerant as the donor parent, such as BRRI Dhan 47, have been released [39]. Although excellent results were obtained by introgressing the Saltol QTL into mega-varieties [36,[39][40][41], this QTL alone is not sufficient to obtain a high degree of salt tolerance under all conditions and stages. A pyramiding approach was proposed several decades ago [22] and Ahmadi et al. [30] showed that even for tolerance at the vegetative stage, the crossing design might involve a large number of accessions that carry different tolerance mechanisms. The question as to which QTLs should be combined must be resolved within each genetic background.
In Europe, because of the risk of low temperatures during the cropping season and the need for accessions to flower during the long days of European summers, almost all rice varieties belong to the temperate japonica group [42]. Thus far, few studies have focused on these accessions, although some studies have evaluated temperate japonicas from other areas [12,29,43]. An association study targeted to specific regions of the genome that carry known tolerance genes was performed [30] but no genome-wide association studies (GWASs) have been conducted thus far. The two objectives of the present study are to evaluate the salinity tolerance at the seedling stage of a panel of 240 temperate japonica rice accessions jointly established by breeders from France, Italy and Spain, and to identify potential targets and candidate genes for marker-assisted breeding using genome-wide association mapping.

Plant materials
The panel comprised 240 temperate japonica accessions. The list of all accessions and their countries of origin is provided in S1 Table. The majority of these lines belong to the panel described in Biscarini et al. [44]. The seeds of these accessions were obtained from the Rice Research Unit of the Consiglio per la ricerca in agricoltura e l'analisi dell'economia agraria (Vercelli, Italy). Seeds from additional Spanish and French lines were obtained from the Instituto de Investigación y Tecnología Agroalimentaria (Torre Marimon, Spain) and from the Centre Français du Riz (Arles, France), respectively. Several accessions known for their susceptibility or resistance to salinity were added to the panel. The indica accessions Nona Bokra, Pokkali, BRRI Dhan 47, FL478, IR64-Saltol, and Hasawi were used as tolerant checks. The recombinant inbred lines BRRI Dhan 47 and FL478 derived from the cross IR29 x Pokkali are recognized for their high salinity tolerance [36]. IR64-Saltol is a version of IR64 developed by the International Rice Research Institute (IRRI) in which the Saltol QTL from Pokkali on chromosome 1 was introgressed by marker-aided selection. IR29 (indica), Giano and Aychade (temperate japonica) were used as susceptible checks. All indica lines were obtained from IRRI.

Genotyping of the panel
The genotyping-by-sequencing (GBS) data used in the present study were derived from pooling the sequences from two batches of accessions (batch 1 and batch 2) obtained from separate libraries. DNA extraction for batch 1 and batch 2 was performed as described in Biscarini et al. [44] and Courtois et al. [45], respectively. Library production and sequencing were conducted in 2013 for batch 1 and in 2016 for batch 2; however, in both cases, the procedures were performed at the Genomic Diversity Facility of Cornell University according to the method of Elshire et al. [46]. The libraries were built using ApeKI for genome digestion. The libraries were sequenced with an Illumina Genome Analyzer II (San Diego, California, USA). The fastq sequences from the two batches were pooled together and aligned to the rice reference genome (Os-Nipponbare-Reference-IRGSP-1.0 [47]). Single nucleotide polymorphism (SNP) calling was performed using the Tassel GBS pipeline v5.2.29 [48]. Five accessions with a rate of missing data above 30% were discarded. Subsequently, markers with a call rate below 80%, a heterozygosity rate above 10% or a minimal number of variant accessions below 5 (corresponding to a minor allele frequency (MAF) < 2.0% in this panel) were discarded. The remaining heterozygotes were converted to missing data. The missing data were imputed using Beagle v4.0 [49]. The final resulting matrix involves 235 individuals and 30 314 markers and can be downloaded in HapMap format from http://tropgenedb.cirad.fr/tropgene/JSP/interface.jsp? module=RICE study "Genotypes", study type "GreenRice_panel_GBS_data".

Phenotyping the panel for salinity tolerance
Experimental design. The experimental design of the trial was a split-plot with three replicates staggered over time. The whole plot main factor was salinity treatment (control (CTRL) versus salt (SALT)) while the secondary split-plot factor was genotype. In each replicate, the plants were distributed into 12 tanks (6 control and 6 salted). The tanks were considered sub-plots. Two resistant checks (Nona Bokra and Pokkali) and three susceptible checks (IR29, Aychade and Giano) were replicated in each tank.
Plant growth conditions. The experiment was conducted in a growth chamber with 28-24˚C day-night temperatures, an 11-13 h day-night photoperiod and a relative humidity of 65%. Tanks of 80 l capacity were used. The tanks were filled with Hoagland and Arnon solution (50 l per tank) [50], pH 5.5, modified as indicated by Courtois et al. [45]. The solution was continuously and gently stirred-aerated using small aquarium pumps (one per tank). The solution was changed every week. The pH was measured every day and adjusted if necessary. Each tank contained a 55.5 cm x 36.6 cm x 3.0 cm foam polystyrene plate with 48 holes drilled at a diameter of 2.5 cm. One hole was left empty and was used for the manipulation of the probe for pH control. Four seeds per accession were sown in the holes of the polystyrene plate. For the first week, a mosquito net with a small mesh was fixed to the bottom of the plate to maintain the seeds. One week after sowing, the polystyrene plates were changed. The most developed seedling of each hole was wrapped in a small 2.5 cm x 2.5 cm x 2.5 cm polyethylene terephthalate fiber cube (Hail mini cubes, Sure to Grow, Ohio, USA) that was cut diagonally through the center. The seed and its attached roots were set below the cube. The wrapped plant was transferred to a new plate without a plastic grid (one plant per hole). Salt treatment was initiated at two weeks after sowing. To limit the osmotic stress, half of the target salt quantity was added to the stressed tanks on day 14. The other half of the target salt quantity was added on day 15. The final concentration was 50 mM NaCl (3 g/l), corresponding to an electrical conductivity of 6.0 dS m -1 . The salt concentration was selected based on the results of Ahmadi et al. [30] to limit lethality while strongly decreasing the growth of the most susceptible accessions. The plants were grown for two weeks under these conditions. The conductivity was adjusted as needed by adding osmosed water. The mortality percentage was assessed at each change of the solution. Four weeks after sowing, the plants were harvested.
Traits evaluated. The effect of the stress on the plant growth rate was measured (Table 1). For each plant, the number of tillers (TIL) was counted. The lengths of the longest leaf (LL) and the longest root (RL) were measured. The blade of the last fully developed leaf was collected, and its length (LGTH) and width (WDTH) were measured. The area of the blade of the last fully developed leaf (LA) was determined as LGTH x WIDTH x 0.75 based on Yoshida's formula [51]. The shoots and roots of each plant were separated and the shoots, roots and last fully developed leaves were oven dried at 72˚C for 4 days, after which the dry matter was weighed (SHOOT, ROOT and LEAF). The specific leaf area (SLA) was calculated by dividing LA by LEAF. The root-to-shoot ratio (R/S) was calculated by dividing ROOT by SHOOT. The whole shoot dry mass of the stressed set of accessions was cut into small pieces and subsequently ground to a powder using a mechanical grinder (one 5 mm-diameter bead per tube; 2 times for 1 min at 1400 rpm). Then, the Na and K mass fractions (Na and K), expressed as a percent of dry matter, were measured for each sample by atomic emission spectroscopy (ICP-AES). The analyses were conducted at the UR59 laboratory in Cirad Montpellier, France (ISO 9001). The Na/K ratio was calculated by division. In addition, several check plants from the unsalted treatment (3 plants of Nona Bokra and 3 plants of IR29 for each replicate) were treated in the same manner. The data are available in S1 Table.

Statistical analyses
Analysis of variance was conducted on the growth parameters using a mixed model with treatment and genotype as fixed effects, and replicate and tank as random effects. The significance of genotype, treatment and genotype x treatment interaction effects were tested. The least square (LS) mean values were calculated and adjusted from the tank effect. For each measured growth trait, considering intrinsic differences in growth rate between accessions, indices of response to stress were calculated based on the LS mean values as iTRAIT = (SALT-CTRL) x 100 / CTRL. Such indices can be strongly affected by the poor growth of some control plants. A Grubbs' test (P < 0.05) was therefore conducted to remove the outliers using XLStat [52]. For the ions measured only under stressed conditions, the analysis of variance was conducted without treatment effects. Broad-sense heritabilities at the genotypic mean level were calculated for each condition and trait by running a separate analysis of variance without treatment effects with h 2 = σ 2 G/(σ 2 G+σ 2 e/n), where σ 2 G and σ 2 e are the estimates of genetic and residual variances, respectively, derived from the expected mean squares of the analysis of variance and n is the number of replicates. The heritabilities of the indices were obtained in the same manner after computing the indices for each replicate. Since some of the measured traits were correlated, a principal component analysis (PCA) was conducted on all traits using XLStat.

Genome-wide association study
A GWAS was conducted on the 231 japonica accessions in the panel having both phenotypic and genotypic data using two different methods, termed the "classical method" and the "interaction method". For the classical method, the GWAS analyses were conducted using a mixed linear model (MLM) with control of panel sub-structure (fixed effect) and kinship among accessions (random effect) using Tassel v5 [53]. The coordinates of the accessions in the two first axes of a PCA conducted on the genotyping data were used to control the population structure. The number of axes was set to two because previous results from an analysis of population structure showed that the main differentiation among European accessions was observed between temperate and tropical-derived accessions [42]. The kinship matrix was calculated using the centered identity by state method, which provides a matrix equivalent to the VanRaden matrix. As described by Courtois et al. [45], a comparison was first conducted between models with kinship alone and kinship and structure using the Akaike information criterion (AIC). The AIC was calculated as -2 ln(L) + 2k, where ln is the natural logarithm, L is the maximized value of the likelihood function for the estimate model and k is the number of estimated parameters [54]. This comparison showed that the second model best fit the data (S2 Table and S3 Fig) and this model was adopted for further analyses. To determine significance and ensure the robustness of the detected associations, a sub-sampling procedure was conducted, according to Tian et al. [55] (maize), and Lafarge et al. [56] and Bettembourg et al [57] (rice). A set of 200 accessions (approximately 80% of the panel size) was selected at random without replacement and a GWAS was conducted on this set. The random sub-sampling and GWAS analysis were repeated 100 times. The number of times an association was detected at P < 0.001 in the 100 sub-samples was counted to obtain a sub-sampling posterior probability for each marker. Based on the distribution of the results, a threshold that had a 95% chance not to be overtaken was selected separately for each trait. The significances reported in the result table for a given marker are those from the GWAS analysis conducted on the entire panel. Markers that were significant within an interval of less than 100 kb were considered to belong to the same QTL. This window size was selected because it was intermediate between the average marker distance (approximately 12 kb in the GBS dataset) and the distance at which the linkage disequilibrium (LD) decayed to half its initial level (estimated at 400 kb in the panel on average across all chromosomes but much more variable in short-range segments). Quantile-quantile (QQ) plots and Manhattan plots were obtained from Tassel or drawn using the qqman R package.
For traits measured both under control and salt conditions, the interaction method, described in detail in Al-Tamimi et al. [33], was also used. The interaction MLM model includes the main effects of the marker and treatment (control and salt, respectively) and marker x treatment interaction, using the same PCA coordinates and kinship matrix obtained from Tassel. This approach enables the identification of SNPs significant for the marker x treatment term and therefore, the identification of loci that respond to salinity treatment. The scripts from Al-Tamini et al. [33] required ASReml-R [58], a proprietary package for the R environment [59]. In this second case, sub-sampling was not attempted because of the time necessary for analyses with ASReml. A false discovery rate (FDR) of 0.05 was used for determining the significance for all traits [60].
In both cases, once a first list of significant markers was established, an LD matrix among these markers was calculated using Tassel, and pairs of close markers with an r 2 value above 0.8 and similar MAFs were pooled in the same interval.

Candidate genes
A list of genes shown to play a role in salt tolerance in rice, referred to as salt tolerance genes, with their identifiers in the two rice gene nomenclatures, physical, position, annotation, function relative to salt stress (when known) and digital object identifier (DOI) of the reference paper(s) was established based on an update of previous lists established by Negrao et al. [13], Kumar et al. [34] and the OGRO database (http://qtaro.abr.affrc.go.jp/ogro). The gene positions and annotations were obtained from Michigan State University (http://rice.plantbiology. msu.edu/), release 7. This list of approximately 300 genes (S2 Table) was used to assess whether the association positions were close to known genes. In addition, all genes with an annotated function in a window of +/-100 kb around the significant markers, or in the interval between two positions when its size was above 200 kb, were extracted using OrygenesDB tools (http:// orygenesdb.cirad.fr/), and their annotations were also explored using key words to assess whether these genes had a plausible relationship with salinity tolerance. The gene HKT1;5 (Os01g20160), located in the center of the Saltol QTL, has a recognized importance in salinity tolerance in rice [14]. To assess the diversity of this gene among temperate accessions, a haplotype analysis of the gene was conducted using the sequences of HKT1;5 from accessions of the 3,000 rice genomes project [61]. The analyzed set comprised the 57 European accessions common to both our panel and the 3,000 genomes project (18 classified as "japx", 30 as "temp" and 9 as "trop" in the IRIC database (http://oryzasnp.org/iric-portal/), with the term "GERVEX" included in their IRGC-ID; listed in S1 Table), 536 random accessions from the 3,000 genomes project, and seven accessions mentioned by Platten et al. [14] as references for the various haplotypes (IR64, Pusa basmati, Pokkali, Nona Bokra, FL478, Nipponbare and Azucena) with known salinity tolerance levels. For these 600 accessions, the sequences involving 76 polymorphic SNPs were extracted from the IRIC database. An unweighted neighbor-joining tree was built using a simple matching index.

Response to salinity
After 28 days of growth with two weeks of salinity stress at a conductivity of 6.0 dS.m-1, the hydroponic experiment was harvested. Lethality affected 4.8% of the control plants and 8.0% of the stressed plants with similar proportions in all replicates, indicating a slight increase in mortality due to salt treatment. However, approximately two-thirds of the plants that died were recorded as having poor germination or were scored as being weak prior to treatment (63.4% for the control plants and 64.7% for the stressed plants).
The ANOVA results showed highly significant treatment and genotype effects and highly significant genotype x treatment interactions for all growth traits except TIL ( Table 2). The significant interactions indicated that the accessions responded differently to the salinity treatment and that variability in tolerance levels existed among the tested accessions. On average, the salinity treatment negatively affected all growth traits, except RL and R/S, which were positively affected (Table 3; Fig 1). LA (-51.2%) and SHOOT (-36.5%) were the traits most impacted by salinity stress, followed by LL (-27.0%) and ROOT (-25.7%). Significant differences between accessions were also observed for Na + and K + leaf mass fractions and the Na/K ratio ( Table 2). Based on the results of the two extreme checks (Nona Bokra and IR29) in which the ion mass fractions were measured under both control and stressed conditions, the mean increase in Na + leaf content varied from 400% (Nona Bokra) to 2000% (IR29), while the mean decrease in K + leaf content varied only from 35% (Nona Bokra) to 56% (IR29). As shown in Fig 2, which is based on the Na/K ratio, the temperate accessions were regularly distributed over an increasing axis from Nona Bokra, Pokkali and BRRI Dhan 47, the most resistant checks with the lowest Na/K ratio (in orange in Fig 2), to Aychade, Giano and IR29, the most susceptible checks with the highest Na/K ratio (in green in Fig 2). No accession was more tolerant than Nona Bokra under the selected stress conditions, but certain accessions, such as Victoria or Rotundus, were close. A few accessions, such as Santerno or YRM6-2, were more susceptible than IR29. When the checks were excluded, the range of variation in the panel was from 0.240% to 1.387% (6-fold difference) for Na, 1.461% to 3.111% for K (2-fold difference) and 0.095 to 0.819 (9-fold difference) for Na/K. The heritabilities of the different traits were moderate to high (0.49 to 0.92) when calculated for each treatment, and for a given trait, the heritabilities were of the same magnitude for both treatments (Table 4). However, the response indices registered much lower heritabilities (0.12 to 0.55). This characteristic translated to the lowest reliability of indices due to propagation of uncertainty in functions involving several variables, and this effect occurred even after removing potential outliers with abnormal growth patterns using a Grubbs' test.
The correlations among indices of growth response to stress were generally high and positive (0.24 to 0.91), except for correlations with iSLA (all negative) and iR/S. The correlation between Na and K was intermediate and negative (-0.39). The correlations between indices of growth response and the Na/K ratio were low and negative; these correlations were significant (-0.23< r <-0.14) for iLL, iROOT, iSHOOT, and iLA, and not significant for iTIL, iR/S and iSLA (S1 Fig), indicating that the phenomena involved in growth response and ion accumulation may be somewhat different. To assess whether the dilution factor played a role in the salinity tolerance of the panel, the correlations between ion-related traits and shoot biomass of the control treatment were calculated; however, the correlations were close to 0 and were not significant. These trends were confirmed by the results of a PCA conducted on the 11 variables. The first two axes explained 57.6% of the variation. The indices of stress response were the main contributors to axis 1, while the ion mass fractions were the main contributors to axis 2 (S2 Fig).

Genome-wide associations
The indica checks were excluded from further analyses, and GWAS analyses were conducted on the 231 accessions that had genotypic data. The eight stress response indices and the three ion parameters were analyzed. The results of the analyses are presented in Table 5 for the classical method and Table 6 for the interaction method. In this narrow genetic base-panel, LD slowly decays. To assess whether the markers, consecutive or not, corresponded to different QTLs, the LD between significant markers was calculated for both methods (S3 Table). No long range or cross-chromosome correlations between significant markers were observed. However, in a few cases, close markers with similar MAFs were found in strong LD (r 2 above For the classical method, a total of 27 independent significant associations validated by subsampling were identified, with the number per trait varying from 1 to 6. Fewer associations were generally detected for the indices of stress response than for the ion concentrations, which is likely due to the lower heritability of the indices than that of the directly measured parameters. The significance was generally moderate (P > 1.0E-06). One association was common between four traits (iLL, iRL, Na and Na/K), and two others associations were common between two traits (Na and Na/K and K and Na/K). The Manhattan plots corresponding to the ions traits are presented in Fig 3 and those for the indices of response are presented in S4 Fig. For the interaction method (Fig 4), a total of 50 significant associations corresponding to qvalues below 0.05 were identified with the number per trait varying from 0 (no QTL detected for LL and R/S) to 15; this range is larger than the range described for the classical method. The significance was also much higher (up to P = 1E-09). More associations were common between traits (four between two traits and five between three traits). These common associations were primarily found between highly correlated traits. Three associations (q01_01, q06_02, and q06_03,) were detected using both methods, among which two associations were detected for the same traits in both methods (q06_02 and q06_03).
The positions of all significant markers were compared to those of the salt tolerance genes listed in S2 Table (Fig 5 for chromosomes 1 to 6 and Fig 6 for chromosomes 7 to 12). A salt tolerance gene was found in the vicinity (< 100 kb) of 23 of the significant markers (13 markers for the classical method and 10 markers for the interaction method). The list includes several ion transporters, notably a sodium transporter (OsHKT1;2), and several genes involved in Ca 2+ signaling or metabolism, including several protein kinases (Table 7). Based on a key word search, the exploration of 1676 genes with tentative functions found in the same intervals (S4 Table) enabled the identification of 84 additional candidates (highlighted in blue in S3 Table). These candidates included 14 antioxidants, 31 transporters, 15 phosphatases, 5 calcium/calmodulin binding proteins, 3 LEA and 3 aquaporins. The remaining candidates were the sole representatives in their broad functional classes.
No associations near the Saltol QTL, located between 10.8 and 11.8 Mb on chromosome 1 and centered on the HKT1;5 gene [36], were significant. The analysis of HKT1;5 polymorphisms of the 57 European temperate accessions in the present panel that were also included in the 3,000 genomes project database showed that the temperate accessions carried haplotypes belonging to closely related branches of the tree. All haplotypes belonged to clusters in which the reference accession did not carry tolerance alleles for this gene (S5 Fig). It is therefore likely that there is no significant association in this zone because no temperate accession carries HKT1;5 tolerant haplotypes.

Discussion
We conducted hydroponic experiments to determine the response of a panel of temperate japonica accessions to mild salinity stress. A mild level of stress (6 dS m -1 ) was selected because of its relevance to observed field conditions in Europe. This mild stress decreased plant growth without inducing lethality, consistent with the results of Ahmadi et al. [30]. The results showed the expected range of variation among the checks; the well-known tolerant lines had the lowest Na/K ratios and the susceptible checks had the highest ratios. The diversity of the panel was evenly distributed between the extreme checks, indicating that some European accessions carried alleles of salinity tolerance that were efficient in situations of mild stress. Although only a mild stress and a single late time point were used, this ranking is useful information for breeders. As shown by Campbell et al. [31] and Al-Tamimi et al. [33], high-throughput image-based phenomics platforms enabling longitudinal assessment might offer a more comprehensive view of the evolution of the response of rice to salt stress over time in further experiments. Chr: chromosome; Pos1-Pos2: limits (bp) of the interval with a LD>0.8; Size: size of the interval in kb; MA: number of accessions carrying the minor allele; iTIL: relative number of tillers; iLL: relative maximum leaf length; iRL: relative maximum root length, iROOT: relative root biomass; iSHOOT: relative shoot biomass; iLA: relative leaf area; iSLA: relative specific leaf area; and iR/S: relative root-to-shoot ratio.
https://doi.org/10.1371/journal.pone.0190964.t005 The GWAS enabled the detection of a set of distinct loci significantly associated with either salinity response indices or Na + or K + mass fractions. Based on our experience in comparing different GWAS approaches in previous datasets, GWASs appear sensitive to small variations in the model (e.g., differences in the number of PCA axes to control the population structure) or in the proposed analysis method for improving the computation speed (e.g., the exact method versus a method that does not involve re-estimating variances for each marker). Small methodological changes translate into different significant associations, although certain associations are generally common between analyses, as observed in maize [62]. The best compromise between computational speed, statistical power and the number of false positives is often challenging to determine for data that are not simulated. Therefore, to best explore the present dataset, we used two different approaches to GWAS analyses, a classical mixed model controlling population structure and kinship [63], and an interaction model advocated by Al-Tamimi et al. [33] for situations, such as in the present study, in which two treatments (stress and control) are compared. For the classical method, which is much less computationally intensive than the interaction method, we attempted a sub-sampling method to explore the robustness of the associations as performed by Tian et al. [55] on maize and Lafarge et al. [56] on rice. Sub-sampling was preferred over bootstrapping because of the existence of a panel structure and the simplicity of file preparation, as a given accession was not repeated several times. Sub-sampling was used to assess the effect of minor variations in MAF in the different samples and, to a lesser extent, examine the effect of population structure. Subsampling enabled us to discard associations that were detected in only in a few samples and occasionally detected in a single sample and had a high probability of being artifacts. Although sub-sampling is a satisfactory method to test the robustness of the detected associations, the time necessary for such computations and subsequent exploitation of the results did not allow the use of more than 100 sub-samples. With the development of appropriate tools, a higher number of samples (500 to 1000) might become manageable, leading to even stronger confidence in the sub-sampling outcome.
We compared the results of the classical method to those of the interaction method that was designed to manage split-plot experiments. The interaction method identified associations with much lower P values than the classical method and a higher number of markers generally supported each peak, as also observed by Al-Tamini et al. [33]. However, although the interaction model involved the control of population structure and kinship using the same matrices as the classical model, the control of false positives was not as good with the interaction method as that observed with the classical method, indicating that part of the increased power is likely obtained at the cost of an inflated rate of false discovery. In addition, the QTLs identified by the interaction method were generally different from those identified using the classical method. Only three QTLs were in common, and these QTLs were not having the lowest P values. This lack of consistency across methods may be partly due to the architecture of the traits Chr  involving a large number of QTLs with small effects, which are borderline in terms of significance. Bian et al. [62] suggested this explanation for the poor correspondence among GWAS results in maize. The lack of consistency can also result from a lack of power in the present study design, although the panel size (230 individuals) is reasonably large for plants and the phenotypic variation within the panel is extensive.
Variations in the associations detected by different methods in the discovery phase imply a need to validate the detected associations, e.g., by replicating experiments with independent panels. However, in addition to the unavoidable GxE interactions in replicating any phenotypic experiment, such replications with independent samples might be particularly difficult in a GWAS context because of the differences in structure, LD, MAF values and marker density between panels [64][65]. These differences can explain why it is challenging to obtain similar results in different GWAS. To strengthen the original evidence, other genetic methods based on gene-based haplotype analyses followed by genetic transformation [66] or on the development of bi-parental mapping populations from panel accessions with complementary haplotypes at significant clusters of markers as proposed by Phung et al. [67] might be better validation choices. However, because of the time and effort that these approaches require, they should be attempted only for the most promising cases. Functional genomic studies, notably gene expression studies, have also been used to provide evidence corroborating or invalidating the results of GWASs [68][69].
Some of the identified associations correspond to or are close to genes involved in salinity tolerance in rice (Table 7). More genes were involved in the response to ionic stress rather than osmotic stress, likely because the measured traits involve only late plant responses. Two genes (OsCBL7 and OsCBL8) belonged to the ionic stress signaling pathway involving calcium. The signaling of ionic stress acts through the Ca 2+ /phospholipase C, salt overly sensitive and calmodulin pathways [34,70]. Genes related to signaling and the signal transduction of the same families were also detected by Al-Tamimi et al. [33] as an early stress response in a GWAS conducted on a panel with different genetic backgrounds (aus and indica). Salt stress induces cytosolic Ca 2+ spiking, which in turn activates Ca 2+ binding proteins, such as calmodulins (e.g., OsCPK1 in Table 7). Two other kinases (e.g., OsSAPK1 and Os07g44330) could also play a role in ionic stress responses.
Seven genes in Table 7 are ion transporters or transcription factors (TFs) that directly activate such genes. HKT1;2 on chromosome 4 (a pseudogene in Nipponbare [71]), is a Na + /K + symporter belonging to the HKT transporter family; the members of this family protect leaves from Na + accumulation by removing Na + from the xylem sap [35]. Several genes were vacuolar antiporter exchange proteins, such as high-capacity cation/H + exchangers (OsCAX1 and OsCAX2), a Mg 2+ /H + exchanger (OsMXH2) or a high-affinity Ca 2+ -ATPase (OsACA6), that play a role in intracellular Ca 2+ -mediated sodium ion homeostasis [72][73]. These genes indicate a preferential mechanism of salt tolerance through Ca 2+ selective accumulation in a temperate background. Zeng et al. [74] demonstrated that Na + -Ca 2+ selectivity was highly correlated with yield under stress and suggested using this parameter as a selection criterion in screening. It may be worthwhile to quantify the Ca mass fraction in future studies involving temperate accessions. OsKAT1 is a highly selective inward-rectifying K + channel protein of the Shaker family, which also plays a role in cellular homeostasis [75]. The OsKAT1 gene was present in a highly significant cluster of SNPs detected in another GWAS conducted in a temperate japonica sub-panel [31]. One gene, OsPCF2, is a TF that binds to the promoter of OsNHX1, which is itself a vacuolar K + -Na + /H + antiporter that plays role in Na + compartmentalization into vacuoles [70].
Four genes (OsNAC1, OsNTL6 (or ONAC070), OsNAC3, and OsNAC5) belong to the large family of NAC (NAM, ATAF and CUC) or NAC-like TFs. Most of these genes are induced by salt stress and regulate salt-responsive genes through an ABA-dependent manner [76]. The role of OsNTL6 and other membrane-bound TFs has been clarified [77]. Environmental changes trigger the release of the TF from the membrane to which it is anchored and, after translocating to the nucleus, this TF regulates the expression of genes involved in stress signaling and stress responses. At least one gene is involved in the scavenging of oxygen reactive species (OsRIP18 with superoxide dismutase activity). OsLEA5 is a late embryogenesis abundant protein. The LEA genes accumulate during water stress, playing a protective role in cellular structures [17]. Two genes belong to the BURP family. Although BURP genes respond to salt stresses, their precise function is unknown [13]. GF14E, OsABF2, and OsHsfB2b, which are isolated members of other families, either regulate or are regulated by salt stress but there is too little information to formulate hypotheses on the precise role of these genes. Among these genes, OsABF2 was found near a QTL for survival days of seedlings in another GWAS with a japonica panel [12].
Although the list in S2 Table includes more than 300 genes, this list is certainly not exhaustive. Therefore, we also examined genes underlying the significant associations/intervals using keywords (S4 Table). Additional potentially interesting genes were observed, but functional evidence of their involvement in salinity response is lacking. Among the transporters, Os06g36590 (a monovalent cation/proton antiporter-2 family), Os03g01330 (a sodium/potassium/calcium exchanger-1), Os07g32530 (a potassium transporter) and Os01g50860 (a chloride transporter), may be relevant. Several genes with antioxidant capacities were identified, and these genes could play a role in the protection of cellular components [13]. Several protein phosphatases, including a large cluster of serine/threonine phosphatases 2C were on the list. Protein phosphatases constitute a large family in rice, and are known to reverse the activity of protein kinases and play a functional role in ABA-mediated signaling pathways and stress responses [78]. Moreover, more than 70 receptor kinases and 80 other kinases were also on the list. Three aquaporins, OsNIP14;1, OsPIP2;7 and OsPIP1;3, were recorded. Aquaporins mediate water transport across cellular membranes and play an important role in maintaining water homeostasis during salinity/osmotic stress [17]. Three LEA proteins were found in addition to OsLEA5 (see above).
The search for candidate genes was fruitful but should not be overvalued. The first limitation of this approach is that, in this panel, full LD (r 2 = 1.0) is often observed between remote points on the chromosomes (e.g., q02_2 or q09_01); therefore, the number of genes underlying these intervals is large, which increases the chances of finding candidates among the salt tolerance genes. Second, regarding the larger list of genes (S4 Table), although some of the candidate genes may be interesting, the functions attributed to these genes are often sufficiently vague that a plausible ex-post rationalization for their association with a trait is generally likely, which is particularly true for TFs and kinases because of their high frequency in the rice genome and large range of putative functions.
The Saltol QTL and the HKT1;5 gene that it carries affect the leaf Na + /K + ratio at the seedling stage while they are not among the segments/genes identified as significant in the present panel. The use of the 3,000 genomes project and the data from Platten et al. [14] enabled us to show that it is likely that none or only a few of the temperate accessions carry a favorable allele at the QTL. Since the panel accessions show a large range of salinity tolerance, it is therefore likely that different mechanisms/genes were selected among the temperate accessions. This result also indicates that some progress in the salinity tolerance of the temperate accessions might be obtained by introgressing Saltol into the temperate background, as attempted in the Neurice project (http://neurice.eu/).

Conclusion
The present study provides information on the salinity tolerance of commonly used accessions in European breeding programs. Through comparisons between GWAS models, the present study provides an idea of the variations in the associations detected by different models and highlights the interest of re-sampling methods in this context. These results emphasize the need to validate putative QTLs and candidate genes by other approaches. However, because of the importance of the associations linked to calcium-mediated ion homeostasis-related genes, the present study also gives interesting clues on which pathways are a priority to explore to obtain a better understanding of the mechanisms of salinity tolerance of temperate rice. iTIL: relative number of tillers; iLL: relative maximum leaf length; iRL: relative maximum root length, iROOT: relative root dry weight; iSHOOT: relative shoot dry weight; iR/S: relative root-to-shoot ratio; iLA: relative leaf area; and iSLA: relative specific leaf area.