Guogen Zhang,Zhiyuan Bi,Jing Jiang,Jingbing Lu,Keyang Li,Di Bai,Xinchen Wang,Xueyu Zhao,Min Li,Xiuqin Zhao,Wensheng Wang,Jianlong Xu,Zhikang Li,Fan Zhang,Yingyao Shi
Genome-wide association and epistasis studies reveal the genetic basis of saline-alkali tolerance at the germination stage in rice
germination stage,genome-wide association study,epistasis,rice,saline-alkali tolerance
Introduction Saline-alkali stress is one of the main abiotic factors limiting rice production worldwide. With the widespread use of rice direct seeding technology, it has become increasingly important to improve rice saline-alkali tolerance at the germination stage. Methods To understand the genetic basis of saline-alkali tolerance and facilitate breeding efforts for developing saline-alkali tolerant rice varieties, the genetic basis of rice saline-alkali tolerance was dissected by phenotyping seven germination-related traits of 736 diverse rice accessions under the saline-alkali stress and control conditions using genome-wide association and epistasis analysis (GWAES). Results Totally, 165 main-effect quantitative trait nucleotides (QTNs) and 124 additional epistatic QTNs were identified as significantly associated with saline-alkali tolerance, which explained a significant portion of the total phenotypic variation of the saline-alkali tolerance traits in the 736 rice accessions. Most of these QTNs were located in genomic regions either harboring saline-alkali tolerance QTNs or known genes for saline-alkali tolerance reported previously. Epistasis as an important genetic basis of rice saline-alkali tolerance was validated by genomic best linear unbiased prediction in which inclusion of both main-effect and epistatic QTNs showed a consistently better prediction accuracy than either main-effect or epistatic QTNs alone. Candidate genes for two pairs of important epistatic QTNs were suggested based on combined evidence from the high-resolution mapping plus their reported molecular functions. The first pair included a glycosyltransferase gene LOC_Os02g51900 (UGT85E1) and an E3 ligase gene LOC_Os04g01490 (OsSIRP4), while the second pair comprised an ethylene-responsive transcriptional factor, AP59 (LOC_Os02g43790), and a Bcl-2-associated athanogene gene, OsBAG1 (LOC_Os09g35630) for salt tolerance. Detailed haplotype analyses at both gene promoter and CDS regions of these candidate genes for important QTNs identified favorable haplotype combinations with large effects on saline-alkali tolerance, which can be used to improve rice saline-alkali tolerance by selective introgression. Discussion Our findings provided saline-alkali tolerant germplasm resources and valuable genetic information to be used in future functional genomic and breeding efforts of rice saline-alkali tolerance at the germination stage.
Genome-wide association and epistasis studies reveal the genetic basis of saline-alkali tolerance at the germination stage in rice

Saline-alkali stress is one of the main abiotic factors limiting rice production worldwide. With the widespread use of rice direct seeding technology, it has become increasingly important to improve rice saline-alkali tolerance at the germination stage.

To understand the genetic basis of saline-alkali tolerance and facilitate breeding efforts for developing saline-alkali tolerant rice varieties, the genetic basis of rice saline-alkali tolerance was dissected by phenotyping seven germination-related traits of 736 diverse rice accessions under the saline-alkali stress and control conditions using genome-wide association and epistasis analysis (GWAES).

Totally, 165 main-effect quantitative trait nucleotides (QTNs) and 124 additional epistatic QTNs were identified as significantly associated with saline-alkali tolerance, which explained a significant portion of the total phenotypic variation of the saline-alkali tolerance traits in the 736 rice accessions. Most of these QTNs were located in genomic regions either harboring saline-alkali tolerance QTNs or known genes for saline-alkali tolerance reported previously. Epistasis as an important genetic basis of rice saline-alkali tolerance was validated by genomic best linear unbiased prediction in which inclusion of both main-effect and epistatic QTNs showed a consistently better prediction accuracy than either main-effect or epistatic QTNs alone. Candidate genes for two pairs of important epistatic QTNs were suggested based on combined evidence from the high-resolution mapping plus their reported molecular functions. The first pair included a glycosyltransferase gene LOC_Os02g51900 (UGT85E1) and an E3 ligase gene LOC_Os04g01490 (OsSIRP4), while the second pair comprised an ethylene-responsive transcriptional factor, AP59 (LOC_Os02g43790), and a Bcl-2-associated athanogene gene, OsBAG1 (LOC_Os09g35630) for salt tolerance. Detailed haplotype analyses at both gene promoter and CDS regions of these candidate genes for important QTNs identified favorable haplotype combinations with large effects on saline-alkali tolerance, which can be used to improve rice saline-alkali tolerance by selective introgression.

Our findings provided saline-alkali tolerant germplasm resources and valuable genetic information to be used in future functional genomic and breeding efforts of rice saline-alkali tolerance at the germination stage.

Soil saline-alkalization, a major factor limiting crop production, has become an increasingly serious global issue in recent years (Shabala, 2013; Wang et al., 2018a). Saline-alkalized soils have two coexisting abiotic stresses, salt stress and alkali stress. Salt stress is caused by neutral salts such as NaCl and Na2SO4, which can mainly cause Na+ toxicity, osmotic stress, and oxidative stress (Zhao et al., 2020). Alkali stress is primarily induced by alkaline salts such as NaHCO3 and Na2CO3 (Fang et al., 2021), resulting in an inhibition of plant growth through additional high pH stress, apart from the effects caused by salt stress (Zhang et al., 2017; Kaiwen et al., 2020). In fact, many reports indicate that alkali stress has more severe effects on crop plants than salt stress (Campbell and Nishio, 2000; Shi and Wang, 2005; Gao et al., 2008). Rice (Oryza sativa L.) is one of the most important food crops and is the primary food source for over half of the world’s population (Huang et al., 2013). Rice is sensitive to saline-alkali stress at different growth stages (Munns and Tester, 2008; Rao et al., 2013). Saline-alkali stress severely disrupts the growth and development of rice plants, resulting in reduced yields and thus is a major threat to world food security (Zhang et al., 2009). Globally, at least 900 million ha of abandoned land are otherwise suitable for crop production because of the severe saline-alkalized soils (Li et al., 2014). In the northeastern and northwestern China, rice growing areas have been expanding rapidly as the global warming. However, rice crops in these areas are facing increasingly more severe saline-alkali problem because of the fast evaporation of irrigated water of rice fields, in addition to ~9 million ha of heavily saline-alkalized wild lands that would have been used for rice production if the saline-alkali problem is overcome. Meanwhile, very limited breeding efforts have been taken to develop and adopt rice varieties tolerant to saline-alkali stress, due partially to “unavailability” of rice germplasm accessions highly tolerant to saline-alkali stresses, and partially to our poor understanding the genetic basis of rice tolerance to saline-alkali stresses. Involving complex physiological and molecular mechanisms (Ganapati et al., 2022), rice saline- alkali tolerance is expected to be controlled by many genes or quantitative trait loci (QTLs). To facilitate breeding efforts for developing saline-alkali tolerant rice varieties for expanding rice growing areas. There are increasing interests to identify the genes and QTLs responsible for saline-alkali tolerance in rice. However, past efforts have been primarily focused on genetic and molecular dissection of salt tolerance (ST) in rice. To date, at least 85 QTLs and genes related to rice salt tolerance have been identified by conventional QTL mapping and map-based cloning from mutants (Yao et al., 2005; Sabouri and Sabouri, 2008; Sabouri et al., 2009; Thomson et al., 2010; Wang et al., 2012; Quan et al., 2017; Dai et al., 2022). According to the funRiceGenes database (Huang et al., 2022), at least 383 genes are associated with salt tolerance, such as SKC1 (Ren et al., 2005), DST (Huang et al., 2009), HST1 (Takagi et al., 2015), qSE3 (He et al., 2019). In contrast, few efforts have been made to understand the genetic basis of alkali tolerance in rice. Li et al. (2020) identified a major QTL, qAT11 on chromosome 11, for alkali tolerance at the bud stage using a population of 184 recombinant inbred lines. Wang et al. (2019) reported that OsPPa6, an osmotic regulatory factor, positively regulates alkaline tolerance rice. Guo et al. (2014) demonstrated that ALT1, encoding the Snf2 family chromatin remodeling ATPase, negatively regulates alkali tolerance mainly by enhancing the defense system against oxidative damage in rice. Genome-wide association studies (GWAS) overcome two main limitations of QTL mapping by traditional biparental genetic populations in plants: (i) providing relatively higher mapping resolution due to plenty of historical recombination events and (ii) using the natural population with more genetic diversity (Korte and Farlow, 2013). Recently, GWAS have become an effective and powerful tool for identifying genes associated with complex traits in rice (Huang et al., 2010; Shi et al., 2017; Liu and Yan, 2019; Zhang et al., 2019; Liu et al., 2021b). In the past decade, GWAS has been successfully used to detect the potential loci for salt tolerance at various developmental stages in rice, which has been well summarized by a recent review (Dai et al., 2022). In contrast, there are few studies on the potential loci for alkali tolerance. Mei et al. (2022) identified 90 loci associated with alkali tolerance at the germination stage by GWAS using 428 diverse rice accessions. Li et al. (2019) detected eight QTLs significantly associated with alkali tolerance-related traits at the seedling stage by GWAS using 295 Geng (japonica) varieties. Although some main-effect loci for saline-alkali tolerance have been identified via GWAS in rice, the genetic architecture of rice saline-alkali tolerance remains unclear since genome-wide epistatic interactions (i.e., QTN by QTN interactions [QQIs]) have been ignored in most GWAS. Epistasis is pervasive and has been recognized as a key aspect of plant genetic architecture (Brachi et al., 2011; Korte and Farlow, 2013; Doust et al., 2014; Mackay, 2014). Recent studies have demonstrated that identifying epistatic interactions in diverse plant germplasm resources can retrieve a large proportion of the “missing heritability”, thereby improving phenotypic prediction accuracy and providing a comprehensive genetic basis for complex traits (Maurer et al., 2015; Luo et al., 2017; Carré et al., 2022). With the rapid population size and marker density increase, GWAS faces statistical and computational challenges in detecting robust epistatic interactions. Fortunately, several new statistical methods and software packages have been developed to identify epistatic interactions in GWAS with high computational efficiency and low false positive rate. Recently, Li et al., 2022a; Li et al., 2022b proposed a powerful and accurate approach, IIIVmrMLM, by integrating the compressed variance component mixed model with the multi-locus random-SNP-effect mixed linear model for detecting QTNs and QQIs in GWAS. To unravel the complex genetic mechanisms underlying saline-alkali tolerance in rice, we conduct a GWAS by the IIIVmrMLM method to detect the genome-wide main-effect QTNs and QQIs associated with seven saline-alkali tolerance-related traits at the germination stage using a large panel of 736 diverse rice accessions from the 3,000 rice genome project (3KRG). By identifying large numbers of QTNs and their QQIs, our findings greatly enriched the current knowledge on the genetic architecture of saline-alkali tolerance in rice and provided valuable germplasm resources for developing saline-alkali tolerant varieties through molecular breeding.

A total of 736 accessions from the 3,000 rice genome project (3KRG) were used to evaluate saline-alkali tolerance of rice at the germination stage. Based on the known population structure (Wang et al., 2018b), these accessions belong to 12 subgroups, including admix (adm) (4 accessions), Aus (160 accessions), Basmati (67 accessions), Geng/japonica (GJ)-adm (8 accessions), GJ-sbtrp (65 accessions), GJ-tmp (58 accessions), GJ-trp (50 accessions), Xian/indica (XI)-adm (78 accessions), XI-1A (56 accessions), XI-1B (59 accessions), XI-2 (61 accessions) and XI-3 (70 accessions) ( Supplementary Table 1 ).

To evaluate saline-alkali tolerance of the rice accessions at the germination stage, we measured seven germination-related traits following the method of Mei et al. (2022). Briefly, surface sterilized seeds that had been broken dormancy were treated with 0.15% Na2CO3 solution (pH = 11) as the saline-alkali stress treatment at the germination stage and with distilled water as the control. Three biological replicates for each accession were performed. In each replicate, 20 seeds of each accession were evenly placed in a 90-mm diameter Petri dish with filter paper, then 10 mL of 0.15% Na2CO3 solution or distilled water was added. We placed the Petri dishes in a growth chamber under the 14-h light/10-h dark photoperiod (28°C/26°C) with 70% relative humidity conditions for 7 d and then changed the Na2CO3 solution and distilled water daily. Germination of seeds was considered to have taken place when the root length was equal to the grain length and the shoot length was equal to half the grain length. We counted germinated seeds under the saline-alkali stress and control conditions every day. The germination potential (GP) and germination rate (GR) of each accession in each replicate were calculated as the proportion of germinated seeds (out of the 20 seeds originally sown) on day 3 and day 7, respectively. Then, germination index (GI) was calculated as GI = Σ(Gt/Tt), and mean germination time (MGT) was calculated as MGT = ΣTtNt/ΣNt, where Gt is the accumulated number of the germinated seeds on day t, Tt was the time corresponding to Gt in days, and Nt was the number of newly germinated seeds on day t (Alvarado et al., 1987; Wang et al., 2010). On day 7, shoot length (SL) and root length (RL) of eight seeds from each accession per replicate in the same growth trend were measured with a ruler. We calculated vigor index (VI) as follows: VI = (mean SL) × GI. The ratios of the saline-alkali stress to the control conditions for the seven germination-related traits, namely relative mean germination time (RMGT), relative germination index (RGI), relative germination potential (RGP), relative germination rate (RGR), relative root length (RRL), relative shoot length (RSL), and relative vigor index (RVI) were used to evaluate rice saline-alkali tolerance during germination.

The 3KRG 1M GWAS SNP dataset was downloaded from the Rice SNP-Seek Database (Alexandrov et al., 2015). We used PLINK 1.9 (Purcell et al., 2007) to obtain a subset of 610,943 SNPs with a missing ratio < 20% and minor allele frequency > 5% for GWAS. The “Single_env” method of IIIVmrMLM (Li et al., 2022a) was used to identify main-effect QTNs associated with saline-alkali tolerance based on a multi-locus GWAS model. The parameters of IIIVmrMLM function were set to “SearchRadius = 20, svpal = 0.01 and num_Threads = 8”. Principal component analysis was performed by GCTA (Yang et al., 2011) with the “-make-grm” parameter and the first five principal components were used to control population structure. The significance thresholds of association by the Bonferroni correction method (0.05/m) were calculated for claiming significant QTNs for the whole panel (P = 8.18 ×10-8), where m is the number of SNPs. The SNP with the minimum P value at a locus was considered the lead SNP. Based on reported genome-wide linkage disequilibrium (LD) decay in 3KRG (Wang et al., 2018b), the other significant SNPs within a range of 150 kb on each side of a lead SNP were merged as a single association locus. The local LD block analysis was performed within 200 kb upstream and downstream of the lead SNP by the “–blocks” function with the default parameter in PLINK 1.9 (Purcell et al., 2007). The heatmap of pairwise LDs was plotted using the LDBlockShow (Dong et al., 2021).

To identify QTN by QTN interactions (QQIs) associated with saline-alkali tolerance-related traits during rice germination, we selected the SNPs with P < 0.01 for each trait except RMGT based on the above multi-locus GWAS. For computational reasons, we screened the SNPs with P < 1×10-5 for RMGT. Then, the “Epistasis” method of IIIVmrMLM was used to detect epistatic QTNs. The parameters of the IIIVmrMLM function were set to “SearchRadius = [0, 1], svpal = [0.10, 0.10] and sblgwas_t = -2.50” and num_Threads = 8”. Principal component analysis was performed by GCTA (Yang et al., 2011) with the “-make-grm” parameter and the first five principal components were used to control population structure. The threshold for significant epistatic QTNs was set at LOD = 3.0. The QQIs results were plotted into circos diagrams using the online website https://www.omicstudio.cn/tool/50 (Gu et al., 2014).

Candidate genes in a QTN locus associated with saline-alkali tolerance were determined based on the following criteria: (1) excluding genes encoding transposon, retrotransposon and hypothetical proteins based on the annotation of Nipponbare reference genome IRGSP 1.0 (Kawahara et al., 2013); (2) significantly regulated (gene expression fold-change > 4.0) by the saline-alkali stress according to Plant Public RNA-seq Database (Yu et al., 2022); and (3) significant phenotypic differences between different gene-haplotypes, which consisted of all SNPs within its 1 kb upstream region and the coding sequence region (ignored synonymous SNPs) (Zhang et al., 2021). Only major haplotypes carried by at least 15 rice accessions were used for phenotypic multiple comparisons by Duncan’s multiple range post-hoc tests with the “agricolae” package in R.

To detect expression levels of the candidate genes of some main-effect and epistatic QTNs, we selected two saline-alkali sensitive accessions (American Huangkedao [Geng] and ARC 6052 [Xian]) and two saline-alkali tolerant accessions (SACHIKAZE [Geng] and ARC 18597 [Xian]) from the GWAS panel. After 24 h of saline-alkali stress with 0.15% Na2CO3, seed embryos of each accession were sampled under saline-alkali stress and control conditions, respectively. Three biological replicates for each accession were performed. Total RNA was extracted from seed embryos using plant RNA extraction kit (Tiangen Biotechnology), and reverse transcribed using reverse transcription kit (Tiangen Biotechnology). Real-time qRT-PCR analyses were conducted with Taq Pro Universal SYBR qPCR Master Mix (Vazyme, Q712-02). UBQ was used as the internal control and the relative expression levels of the target genes were calculated using the 2-ΔΔCT method (Livak and Schmittgen, 2001). All primers used for qRT-PCR are listed in Supplementary Table 2 .

To determine the impact of epistatic QTNs on genomic prediction, we compared the prediction abilities based on different QTN datasets for each trait using the R package “rrBLUP” with the “mixed.solve” function (Endelman, 2011). The prediction abilities of the main-effect QTNs, QQIs and their combinations for each trait were obtained based on 10 times repeated 10-fold cross-validation, respectively. The missing genotypes of QTNs were replaced by their mean values. Two-sided Wilcoxon rank-sum tests were used to test the differences among the prediction abilities.

One-way ANOVA was used to test the phenotypic differences for measured traits between subpopulations by the R package “agricolae”. Spearman’s correlation coefficients between different saline-alkali tolerance traits were calculated using the R package “ggcor”. Boxplots and barplots on the measured trait variation in different rice subpopulations were generated using the R “ggplot2” package. The haplotypic enrichment analysis in a subpopulation was performed using the function “phyper” in R.

To evaluate saline-alkali tolerance of the 736 diverse rice accessions at the germination stage ( Figure 1A ), we measured seven germination-related traits (MGT, GI, GP, GR, RL, SL, and VI) under control and saline-alkali stress conditions. As shown in Figure 1B , GI, GP, GR, RL, SL, and VI were significantly reduced and MGT was significantly prolonged under saline-alkali stress compared to the control, indicating that saline-alkali stress significantly inhibited rice germination and growth and the measured traits were indeed related to rice saline-alkali tolerance. The degree of saline-alkali damage to these germination-related traits was as follows: RL, GP, VI, GI, GR, SL, and MGT ( Figure 1C ). Of these, RL, the trait most affected by saline-alkali stress, was reduced by an average of 73.1% across the 736 accessions. To eliminate the influence of genotype on the development of germination traits, we used the relative trait value (stress/control) to evaluate the saline-alkali tolerance of rice accessions. A wide range of variation in the saline-alkali tolerance traits was observed among different rice populations and among different accessions within each of the populations ( Figure 1D ). At the population level, the overall rank of saline-alkali tolerance was Basmati > Xian > Geng > Aus, which varied slightly depending on specific saline-alkali tolerance traits ( Figure 1D ). There was no significant difference in RMGT among the four subpopulations. For the other six traits, the Basmati accessions had the highest mean values in four traits except for RGP and RSL, and the Aus accessions had the lowest mean values in five traits except for RGP. The Xian accessions had significantly higher RGP and RSL than the Geng accessions ( Figure 1D ). Among the three Geng subgroups, GJ-tmp accessions had significantly higher RGP and RSL than the GJ-sbtrp and GJ-trp accessions, while the GJ-trp accessions had significantly higher RRL than the GJ-sbtrp and GJ-tmp accessions. Among the four Xian subgroups, the Xian-3 subgroup had the strongest saline-alkali tolerance, while Xian-1B had the weakest saline-alkali tolerance ( Supplementary Figure 1 ). These findings suggest that identifying loci associated with saline-alkali tolerance and their corresponding favorable alleles from diverse subpopulations is valuable in molecular breeding for rice saline-alkali tolerance. Of the 736 accessions evaluated, 30 accessions from ten countries showing high levels of saline-alkali tolerance at the germination stage were screened according to the two traits of RRL and RGP, including 7 Basmati accessions, 21 Xian accessions, 1 Aus accession and 1 Geng accession, respectively ( Supplementary Table 1 ), providing a valuable source of materials for future theoretical and breeding research on saline-alkali tolerance of rice. As expected, moderate negative correlation existed between the measured saline-alkali tolerance-related traits (RGR, RGI, RVI, RSL and RRL) and RMGT, while positive correlation existed between RGR, RGI, RVI, RSL and RRL, though much stronger between RGI, RGR and RVI ( Figure 1E ). Similar results were found in different subpopulations, though the trait correlation was stronger in population Xian than in other populations ( Supplementary Figure 2 ).

Using the IIIVmrMLM software, 165 QTNs associated with saline-alkali tolerance-related traits at the germination stage were identified based on a multiple-locus GWAS ( Supplementary Table 3 ). Of these, 32, 39, 29, 38, 22, 31 and 25 QTNs were detected for RMGT, RGI, RGP, RGR, RRL, RSL, and RVI, respectively, which explained 49.0%, 63.5%, 52.1%, 59.0%, 38.7%, 61.0% and 44.8% of total trait phenotypic variation, respectively ( Figures 2A, B ). A total of 70 (42.4%) main-effect QTNs co-localized with previously reported QTLs related to alkali tolerance and genes related to saline-alkali tolerance in rice, suggesting that most QTNs detected by the multi-locus random-SNP-effect mixed linear model were reliable ( Figure 2C ). The large number of main-effect QTNs identified in this study also provided evidence for the complexity of the genetic mechanisms underlying saline-alkali tolerance during rice germination. Interestingly, the same QTNs were significantly associated with two or more saline-alkali tolerance-related traits (RGR, RGI, RVI, RSL and RRL) in 29 cases ( Figure 2D ). For example, a main-effect QTN on chromosome 5, q5.7, was commonly associated with RGR, RGI, RSL, and RRL, explaining a small portion of phenotypic variation (2.69% for RGR, 4.68% for RGI, 3.43% for RSL, and 2.59% for RRL). In addition, the most significant main-effect QTN associated with RGR, q1.11, explained 4.49% of the phenotypic variation of RGR.

According to the criteria for candidate gene prediction (see Methods), we identified four and two candidate genes in the important main-effect QTNs, q1.11 and q5.7, respectively ( Figure 2E ). q1.11 was the most significant (P =5.10 × 10-80) QTN associated with RGR. The LD block region of q1.11 was estimated from 25.820 to 25.866 Mb (46.3 kb interval) on chromosome 1 ( Figure 3A ). The lead SNP associated with RGR, rs1_25863302, is located in the promoter region of a putative auxin efflux carrier gene LOC_Os01g45550 (OsPIN3t). Five major haplotypes (n ≥ 15 accessions) were detected using ten SNPs in the 1 kb upstream region of LOC_Os01g45550 ( Figure 3B ). Among the five haplotypes, Hap4 with the largest mean RGR (0.85) was considered the favorable haplotype for saline-alkali tolerance ( Figure 3D ), which was present only in Geng accessions ( Figure 3E ). OsPIN3t was significantly up-regulated in saline-alkali-sensitive Geng accession (American Huangkedao) and saline-alkali-tolerant Xian accession (ARC 18597) under saline-alkali stress, whereas it was significantly down-regulated in saline-alkali-tolerant Geng accession (SACHIKAZE) under saline-alkali stress ( Figure 3C ). Interestingly, American Huangkedao carries the non-favorable haplotype Hap3 and SACHIKAZE carries the favorable haplotype Hap4. In contrast, q5.7 was the most significant (P = 1.94 × 10-30) one associated with four of five saline-alkali tolerance traits (RGR, RGI, RVI, RSL and RRL) except RVI, while the q3.6 (P = 4.71 × 10-12) and q11.16 (P = 8.12 × 10-13) were the QTNs associated with four saline-alkali tolerance traits except RGI and RSL, respectively. The total phenotypic variance explained by this QTN was 2.69% for RGR, 4.68% for RGI, 3.43% for RSL, and 2.59% for RRL. The q5.7 LD block region of lead SNP (rs5_14379574) associated with RGR was estimated from 14.206 to 14.406 Mb (200.0 kb interval) on chromosome 5 ( Figure 3F ). According to reported information of gene annotations ( Figure 2E ), we found that there was a potential salt-tolerant gene LOC_Os05g24770 encoding a reticulon domain-containing protein in the LD block region. LOC_Os05g24770 was significantly up-regulated in saline-alkali-sensitive Geng accession (American Huangkedao), while significantly down-regulated in saline-alkali-tolerant Geng accession (SACHIKAZE) under saline-alkali stress ( Figure 3H ). Three major haplotypes (n ≥ 15 accessions) of LOC_Os05g24770 were detected based on 13 SNPs in its 1 kb upstream region and the coding sequence region ( Figure 3G ). Among the three haplotypes, Hap2 had the largest mean RGR of 0.82, while the mean RGR of Hap1 and Hap3 were 0.77 and 0.62 ( Figure 3I ). Thus, we defined Hap2 as the favorable haplotype for saline-alkali tolerance, which was enriched in Xian accessions ( Figure 3J ).

Using the IIIVmrMLM software, we identified 134 significant QQI pairs ( Supplementary Table 4 ; Figure 4A ) for the saline-alkali tolerance-related traits between 204 epistatic QTNs ( Supplementary Table 5 ), including 80 (39.2%) of the identified main-effect QTNs ( Figure 4G ). The identified QQI pairs (explained total trait phenotypic variance) were 26 (60.0%) for RGP, 15 (31.1%) for RMGT, 22 (39.5%) for RGI, 26 (43.6%) for RGR, 23 (47.0%) for RRL, 15 (27.7%) for RVI, and 11 (32.1%) for RSL, respectively ( Figures 4B, C ). Most (90.5%) epistatic QTNs were involved in single pairwise epistasis, and were associated with only a single trait ( Figure 4D ). Interestingly, 33 and 10 epistatic QTNs associated with two and three traits, respectively ( Figure 4H ). This result indicated that most epistatic QTNs were highly specific regarding the QQIs and traits involved. Notably, 79 (38.7%) epistatic QTNs co-localized with the previously reported QTLs related to alkali tolerance and genes related to saline-alkali tolerance in rice ( Figure 4F ), while the QQI pair was detected between rs3_17097612 and rs7_8490713 with the largest PVE (11.3%) on RRL. The QQIs detected in this study could be divided into three types: (I) QQIs occurred between two main-effect QTNs; (II) QQIs occurred between one main-effect QTN and one non-main-effect QTN; and (III) QQIs occurred between two non-main-effect QTNs. The proportion of different QQI types varied considerably across different traits ( Figure 4E ). Specifically, the proportion of type I QQIs was highest for RGI (50.0%), followed by RVI (33.3%), RMGT (33.3%), RRL (21.7%), RGR (19.2%), RSL (18.2%) and lowest (11.5%) in RGP. The frequency of type II QQIs was highest (60.0%) for RMGT, followed by RGP (57.7%), RVI (53.3%), RGR (50.0%), RRL (47.8%), RSL (45.5%) and lowest (31.8%) in RGI. In contrast, the proportion of type III QQIs was highest (36.4%) for RSL, followed by RGP (30.8%), RGR (30.8%), RRL (30.4%), RGI (18.2%), RVI (13.3%) and lowest (6.7%) in RMGT ( Figure 4E ).

Two important of QQIs with large PVE were noted. The first case was the QQI between rs2_31646695 and rs4_436064, which explained 4.43% of the phenotypic variation for RGP. The mean phenotypic values ​​of the two alleles (AA and TT) at rs2_31646695 were significantly different, with the mean RGP values ​​of AA and TT being 0.45 and 0.24, respectively ( Figure 5A ), while the phenotypic values of the two alleles (AA and GG) at rs4_436064 were also significantly different, with the mean RGP values ​​of AA and GG being 0.46 and 0.33, respectively ( Figure 5B ). The phenotype of the accessions with the non-favorable allele GG at rs4_436064 was further enhanced by the favorable allele AA at rs2_31646695 (i.e., the mean RGP increased from 0.24 in GGTT to 0.37 in GGAA) ( Figure 5C ). Based on the literature search and expression profiles under abiotic stress conditions, two promising candidate genes of the QQI between rs2_31646695 and rs4_436064 were identified, including an E3 ligase gene LOC_Os04g01490 (OsSIRP4) near rs4_436064 and a glycosyltransferase gene LOC_Os02g51900 (UGT85E1) near rs2_31646695. UGT85E1 was significantly down-regulated in saline-alkali-tolerant Geng accession (SACHIKAZE) and saline-alkali-sensitive Xian accession (ARC 6052) under saline-alkali stress ( Figure 5D ). OsSIRP4 was significantly up-regulated in saline-alkali-sensitive Geng accession (American Huangkedao) and saline-alkali-sensitive Xian accession (ARC 6052), while significantly down-regulated in the saline-alkali-tolerant Xian accession (ARC 18597) under saline-alkali stress ( Figure 5E ). Based on the gene-haplotype analysis using the SNPs in the 1 kb upstream region and the coding sequence region, we found significant differences in RGP between the favorable and non-favorable haplotypes of both UGT85E1 and OsSIRP4 ( Figures 5F, G ). Moreover, two haplotype combinations of the two genes were obtained in the GWAS panel ( Figure 5H ). The favorable combination consisted of two favorable haplotypes of the two genes, with an average RGP value of 0.71 ( Figure 5H ). The favorable haplotype of UGT85E1 was mainly enriched in population Xian, while the favorable haplotype of OsSIRP4 and the favorable haplotype combination were carried only by Xian accessions ( Figures 5I–K ). The second case was the QQI between rs2_26441823 and rs9_20644583, explaining 4.52% of the phenotypic variation of RGP. The phenotypic values ​​of the two alleles (CC and TT) of rs2_26441823 were significantly different, with the mean RGP values ​​of CC and TT being 0.43 and 0.34, respectively ( Figure 6A ), while the phenotypic values of the two alleles (AA and TT) at rs9_20644583 were also significantly different, with the mean RGP values ​​of AA and TT being 0.24 and 0.41, respectively ( Figure 6B ). The phenotype of the accessions with the non-favorable allele TT at rs2_26441823 was further enhanced by the favorable allele TT at rs9_20644583 (i.e., the mean RGP increased from 0.12 in TTAA to 0.37 in TTTT) ( Figure 6C ). Based on the literature search and expression profiles under abiotic stress conditions, two promising candidate genes of the QQI between rs2_26441823 and rs9_20644583 were identified, including an ethylene-responsive transcriptional factor gene LOC_Os02g43790 (AP59/OsBIERF3) near rs2_26441823 and a Bcl-2-associated athanogene protein gene LOC_Os09g35630 (OsBAG1) near rs9_20644583. AP59 was significantly down-regulated in saline-alkali-tolerant Geng accession (SACHIKAZE) and saline-alkali-tolerant Xian accession (ARC 18597) under saline-alkali stress, while significantly up-regulated in saline-alkali-sensitive Xian accession (ARC 6052) under saline-alkali stress ( Figure 6D ). OsBAG1 was significantly up-regulated in saline-alkali-sensitive Geng accession (American Huangkedao) and saline-alkali-tolerant Geng accession (SACHIKAZE) under saline-alkali stress, while significantly down-regulated in the saline-alkali-tolerant Xian accession (ARC 18597) under saline-alkali stress ( Figure 6E ). Significant differences in RGP were both observed between the favorable and non-favorable haplotypes of AP59 and OsBAG1 ( Figures 6F, G ). Among the two haplotype combinations of the two genes in the GWAS panel ( Figure 6H ), the favorable saline-alkali tolerant combination was made up of the two favorable haplotypes from the two genes, with an average RGP value of 0.46 ( Figure 6H ). The favorable haplotypes of AP59 and OsBAG1 were mainly carried by Xian accessions ( Figures 6I, J ), while the favorable haplotype combination was only enriched in population Xian ( Figure 6K ).

To determine whether the identified epistatic QTNs have advantages for genomic prediction of saline-alkali tolerance, we compared the performances of the main-effect QTNs, the epistatic QTNs and their combined datasets in genomic prediction for the seven saline-alkali tolerance-related traits using ridge regression best linear unbiased prediction (rrBLUP) approaches. Based on 10 times repeated 10-fold cross-validations, the main-effect QTN dataset showed significantly higher prediction abilities than the epistatic QTN dataset for all traits except for RGP ( Figure 7 ). Notably, the combined dataset of main-effect and epistatic QTNs showed significantly higher prediction abilities than either main-effect or epistatic QTNs alone for all traits. In other words, the prediction abilities of the seven traits based on the main-effect QTNs were significantly improved by inclusion of the epistatic QTNs, with an average increased power of 7.0%, when the epistatic QTNs were included in the prediction model. Among the seven traits, the prediction ability of RGP was improved most (with an increase of 14.7%) by the epistatic QTNs compared to the prediction ability by the main-effect QTNs, while the prediction ability of RMGT was improved least (with an increase of 1.1%) ( Figure 7 ).

As the global climate changes, soil saline-alkalization has become an increasingly serious threat to world crop production and food security. To facilitate the development of rice cultivars with improved saline-alkali tolerance, we aimed to dissect the genetic architecture of rice saline-alkali tolerance using a modified strategy of GWAS. The identification of 165 main-effect QTNs and 134 QQIs that explained an average of 52.6% and 40.1% total PVE for each of the rice saline-alkali tolerance traits during germination was consistent with the current knowledge that rice saline-alkali tolerance is a complex trait involving large numbers of loci and multiple genetic mechanisms (Zhang et al., 2013). Of these main-effect QTNs, 70 (42.4%) overlapped with previously reported QTLs/genes related to saline-alkali tolerance in rice, indicating a very low false positive rate in QTN detection of this study. Moreover, we identified two promising candidate genes, LOC_Os01g45550 (OsPIN3t) and LOC_Os05g24770, associated with saline-alkali tolerance at the germination stage in rice. LOC_Os01g45550 (OsPIN3t) is involved in drought stress response in rice and acts in auxin polar transport (Zhang et al., 2012). Both drought and saline-alkali stresses are known to cause osmotic stress in plant root cells. Also, salt stress is known to significantly impact auxin transport in plants, which is expected to influence plant root development through its effect on auxin distribution and epidermal biosynthesis (Korver et al., 2018). LOC_Os05g24770 encodes a reticulon domain-containing protein in rice. The reticulon domain-containing protein gene, F2DHE3 in Tibetan wild barley was reportedly involved in salt adaptation of the wild barley as it was up-regulated in the variety under high Na+ accumulation (Shen et al., 2016). Furthermore, another reticulon family protein gene, RTNLB2, in Arabidopsis was also up-regulated under salt stress (Ma and Bohnert, 2007). In both cases, OsPIN3t and LOC_Os05g24770 are suggested to be involved in saline-alkali tolerance, but how they are contributing to saline-alkali tolerance remains to be experimentally validated. One major discovery of this study was the importance of epistasis in determining rice saline-alkali tolerance. Although epistasis is known to be an important genetic basis of complex traits (Yu et al., 1997; Li et al., 2001; Marchini et al., 2005; Evans et al., 2006; Wei et al., 2014), it remains a huge challenge to detect epistasis in natural populations by GWAS (Manicacci et al., 2009; Wurschum et al., 2011; Zhang et al., 2015; Wen et al., 2016; Luo et al., 2017). This is because the power in detecting QTNs, particularly epistatic QTNs, is influenced much more by the population structure of the sampled accessions in a GWAS study, which results from the presence of multiple alleles at virtually all loci in a random germplasm population, and varied frequencies of different alleles and multi-locus genotypes at specific loci across different populations of the sample accessions. Fortunately, the newly developed method by Li et al., 2022a; Li et al., 2022b made it possible to detect epistasis by GWAS. For the two cases of important QQIs with large PVE ( Figures 5 , 6 ), we identified four promising candidate genes. For the first case of important QQI between rs2_31646695 and rs4_436064, LOC_Os02g51900 (UGT85E1) and LOC_Os04g01490 (OsSIRP4) were identified as the promising candidate genes. OsSIRP4, highly induced under salt stress in rice, acts as a negative regulator in the plant response to salt stress via the 26S proteasomal system regulation of substrate proteins, OsPEX11-1, which provides important information for the adaptation and regulation in rice (Kim and Jang, 2021). The expression of UGT85E1 was significantly up-regulated by drought stress and ABA treatments, while its overexpression led to the enhanced tolerance in transgenic rice plants to drought stress (Liu et al., 2021a). Salt or drought stress generates osmotic stress in plant cells and dramatically increases cellular ABA levels (Fernando and Schroeder, 2016). For the second case of important QQI between rs2_26441823 and rs9_20644583, LOC_Os02g43790 (AP59/OsBIERF3) and LOC_Os09g35630 (OsBAG1) were identified as the promising candidate genes. AP59/OsBIERF3 was up-regulated by salt, cold, drought and wounding in rice (Cao et al., 2006). The overexpression of AP59 in rice under the control of the constitutive promoter OsCc1 increased the tolerance to drought and high salinity at the vegetative stage (Oh et al., 2009). OsBAG1 was up-regulated (at least four-fold) after continuous high-temperature treatment and was down-regulated in two varieties (Dagad deshi and IR20) under drought stress (Zhou et al., 2021). As heat, drought and salt stress are the main abiotic stresses, we speculate that OsBAG1 may also respond to salt stress. The function of these promising candidate genes of important QQIs should be experimentally validated in future. In fact, this study, to our knowledge, was the first effort to detect genome-wide epistasis for saline-alkali tolerance by GWAS using natural rice germplasm populations. Nevertheless, two interesting observations were noted regarding how to detect epistasis using GWAS of natural populations and the interpretation of the detected epistatic QTNs. Firstly, we detected more epistatic QTNs than main-effect QTNs, and 79 (38.7%) of the identified epistatic QTNs overlapped with previously reported QTLs/genes related to saline-alkali tolerance in rice, indicating a very low false positive rate in the epistatic QTNs and the reliability of the method used in this study, even though most of the identified QQIs appeared to have smaller effects than the main-effect QTNs. Secondly, significantly more (39.2%) main-effect QTNs were involved in QQIs in this study than in previous reports from segregating populations of bi-parental crosses in which most QQIs occur between QTNs that do not have detectable main-effects (Li et al., 1997; Li et al., 2001). Theoretically, this was expected because most QQIs were undetectable in natural germplasm populations due to the fact that the haplotype combination frequencies of any specific pair of QQI in the 3KRG were much lower than their Mendelian segregation frequencies in bi-parental segregating populations ( Figures 5K , 6K ). The impact of the population structure on the power in detecting epistasis is particularly more pronounced in inbreeding species such as rice than in outcrossing ones (Doust et al., 2014). Thus, the epistatic QTNs affecting saline-alkali resistance detected in this study were important ones, which can be verified genetically using bi-parental or multiple-parent mapping populations with balanced frequencies of bi-allelic or multi-locus genotypes in which the epistatic variance is maximized (Mackay, 2014; Sakai et al., 2021; Wang et al., 2022). More importantly, molecular validation of epistatic QTNs identified by GWAS will be facilitated by the higher resolution of their candidate genes, and by their functional relationships of the candidate genes, i.e., to experimentally validate if there is one-way (regulating or regulated) or mutual functional dependency between the interacting QTN candidates (Zhang et al., 2011). Taking together, our results supported that epistasis is an important genetic basis of saline-alkali tolerance in rice.

The ultimate goal of past theoretical work of rice, including both genetic dissection of complex traits and functional genomic research aiming at cloning and functional characterization of genes controlling complex traits, is to apply these results to accelerate genetic gains in improving complex traits in future breeding efforts. One common strategy is to select and pyramid favorable alleles at many target loci to improve traits such as yield potential and abiotic stress tolerances in breeding. It is also clear that virtually all complex traits are controlled by complex signaling pathways involving large numbers of loci. However, the presence of epistasis would suggest that it is unreliable to predict the phenotypic effects of breeding progeny based on their genotypes at single loci. Thus, it is crucial to discover the epistatic interactions between genes to optimize gene pyramiding experiments (Yan et al., 2011). In these cases, taking both the additive and epistatic effects of alleles among multiple genes into consideration would allow identification of the most favorable allele combination(s) at a specific set of epistatic loci with maximum effects on target traits, which otherwise can’t be achieved through simply pyramiding multiple alleles at single loci (Benfey and Mitchell-Olds, 2008). This was well demonstrated in our genome-wide selection experiments in which inclusion of both main-effect and epistatic QTNs had consistently higher predictability for all saline-alkali resistance traits. Similar results were achieved in an effort to improve rice saline-alkali tolerance at the germination stage by pyramiding the favorable haplotypes of multiple saline-alkali-tolerant main-effect QTNs and creating the optimal haplotype combinations (Mei et al., 2022). Obviously, our results suggest that the identified epistatic QQIs and their frequency distribution in different subpopulations can provide valuable information for molecular breeding for improving complex traits including rice saline-alkali tolerance. Thus, it is now possible to enhance saline-alkali tolerance of saline-alkali sensitive rice varieties by introgressing the favorable haplotype combinations via marker-assisted selection and chip-based genotyping. For example, the favorable haplotype combinations, UGT85E1 Hap1 × OsSIRP4 Hap2 and AP59 Hap1 × OsBAG1 Hap1 of the epistatic QQIs can be used to improve RGP by selective introgression ( Supplementary Table 6 ). Clearly, the genetic information of the germplasm resources carrying favorable haplotype combinations obtained in this study would provide useful information for selecting suitable genetic materials for further functional genomic research, though their values for improving rice saline-alkali tolerance should be actively tested in real breeding efforts in future.

Through high-density SNPs, extensive phenotyping of the large-scale rice natural population under the saline-alkali stress and control conditions, and recent advancement in GWAS methodologies, we were able to conduct more robust and efficient genetic analysis of rice saline-alkali tolerance to better address QQI effects. Our results suggest that rice saline-alkali tolerance is a very complex trait with epistasis as an important genetic basis. The epistatic QQIs and their frequency distribution in germplasm resources can provide valuable information for future molecular breeding for improving rice saline-alkali tolerance at the germination stage.

The original contributions presented in the study are included in the article/ Supplementary Material . Further inquiries can be directed to the corresponding authors.

FZ and YS conceived, designed the experiments; GZ, JJ, JL, KL, DB, XW and XYZ performed the experiments; GZ, ZB, ML, XQZ, WW and JX performed analysis and interpretation of the data; GZ, FZ and ZL wrote the manuscript; FZ revised the manuscript. All authors contributed to the article and approved the submitted version.