Rice (Oryza sativa L.) is the most important staple crop of more than half of the global population (Zhao et al. 2017; Pandey et al. 2002). It is widely grown in varied environ-mental conditions, ranging from sea-coast to high altitude. Rice is one of the most flood-threatened crops, with almost 15 million hectares of rice fields in South and Southeast Asia regularly affected by submergence (Neeraja et al. 2007). Oladosu et al. (2020) reported about 700 million people living in poverty in Asian countries who reside in flood-prone rice-cultivating regions of South Asia. The increasing submergence threat associated with climate change and requirements to feed the human population represent a great challenge to increasing potential yields in current agriculture (Godfray et al. 2010). Therefore, the development of crops resistant to submergence is required to solve this problem.

Submergence refers to the condition where the entire aerial part of the plant is completely submerged underwater. This condition can cover the entire plant for prolonged periods over 7 days, which is lethal to most rice varieties. During submergence, plants experience low oxygen con-ditions due to the impeded diffusion of gases and may encounter hypoxia, or in severe cases, even anoxia (Bailey-Serres et al. 2010; Winkel et al. 2013; Mondal et al. 2020; Mital et al. 2022).

To cope evolutionarily with oxygen deficiency in sub-mergence conditions, rice acquired adaptive mechanisms which reflect the ability of the coleoptile to emerge and grow quickly to contact with air (Ismail et al. 2009). According to Stünzi and Kende (1989), submergence resistance varieties were reported to have the ability during submergence conditions to elongate at rates that are 25-fold higher than air-grown plants. Also, in these varieties, if the energy source required for growth is limited due to a lack of oxygen in the submergence condition, it is reported that oxygen is transferred to the roots and drainage channels to survive by rapidly growing the early leaves to obtain the oxygen needed for survival (Hwang et al. 1999).

To minimize the damage caused by submergence, the International Rice Research Institute (IRRI) developed varieties (Swarna-Sub1, IR64-Sub1, Samba Mahsuri-Sub1, BR11-Sub1, TDK1-Sub1, and CR1009-Sub1) with SUB1 that had positive effects for increased submergence tolerance (Septiningsih et al. 2012; Iftekharuddaula et al. 2016). In this manner, varieties to develop breeding lines to enhance tolerance are essential to capture the genetic complexity of the trait and provide quantitative trait loci (QTLs) targets for marker-assisted breeding (Angaji et al. 2010).

Septiningsih et al. (2015) investigated the F_2:3 population between IR72 and Madabaru, both moderately tolerant varieties. In this population, phenotypic data showed trans-gressive segregation with several families having an even higher survival rate; four QTLs were identified on chromo-somes 1, 2, 9, and 12. Gonzaga et al. (2017) conducted QTL analysis to identify QTLs associated with submergence tolerance using 115 F₇ recombinant inbred lines (RILs) derived from the cross of Ciherang-Sub1. The 6K single nucleotide polymorphism (SNP) chip was utilized for QTL mapping, which resulted in 469 polymorphic markers, and two QTLs associated with submergence tolerance were identified on chromosomes 2 and 8. Nghi et al. (2019) conducted a genome-wide association study (GWAS) analysis to detect the QTL region for submergence tolerance. A panel of 273 japonica rice accessions was used for the analysis and discovered 11 significant marker-trait associations and identified candidate genes potentially involved in coleoptile length.

GWAS based on linkage disequilibrium (LD) and single sequence repeat (SSR) markers or SNP markers has been widely performed to identify loci significantly associated with many traits in model plant species, including rice (Ham et al. 2022).

To provide genetic information on submergence toler-ance for rice breeders and geneticists, we here conducted a GWAS to suggest candidate genes that are associated with the submergence tolerance trait.

Plant materials and genomic data

A core set of 137 rice accessions from the National Agrobiodiversity Center of the Rural Development Ad-ministration (RDA, Korea) was used to detect variations in submergence tolerance (Supplementary Table S1). The 137 rice accessions were suggested by Kim et al. (2007) as the Korean Rice Core Selection, which had been collected from 28 countries, selected from 25,604 rice accessions, and could be divided into 6 subgroups: tropical japonica (19 accessions), temperate japonica (62), indica (43), aus (8), aromatic (3), and admixture (2) (Kim et al. 2007; Zhao et al. 2010). The genomic data of the 137 rice accessions were obtained from an average coverage of approximately 8X on an Illumina HiSeq 2500 Sequencing Systems Platform (Illumina Inc., San Diego, USA). Raw reads were aligned against the rice reference genome (IRGSP 1.0) for genotype calling. To generate the genotype dataset, the following parameters were used for GWAS: missing value < 1%, a minor allele frequency (MAF) > 5%, and a heterozygosis ratio < 5%, as implemented using Plink software (Purcell et al. 2007). Finally, approximately 2 million high quality SNPs were obtained from 6.5 million raw data SNPs for the further GWAS (Kawahara et al. 2013; Kim et al. 2016).

Evaluation of survival rate in the submergence condition

After filling the soil in the plastic sowing box (50 holes), 30 dry seeds per variety were divided into three sowing holes by 10 dry seeds, respectively, and were sown to a thickness of 1 cm. The investigation of submergence tolerance was performed using the Kim et al. (2019) and Jang et al. (2022) methods. The plastic sowing box was then moved to a water storage box (length × width × depth; 67 cm × 50.5 cm × 41.5 cm) for water, and the depth of the water was maintained during the test period after water treatment to a depth of 10 cm. The submergence treatment was performed in a growth chamber with an average temperature set to 21℃, and the survival rate was calculated for the individual with leaves appearing on the surface of the water after 21 days of seeding.

Population and genotype analysis

The population structure and cross-validation (CV) analy-sis of the 137 accessions were analyzed using ADMIXTURE version 1.3.0 (Alexander et al. 2009), with subgroups assigned according to the delta K value. For each result, visualizations were performed in the Structure Plot V2.0 web application (http://omicsspeaks.com/strplot2/, accessed on 1 November 2022), and the R calculation of principal component analysis (PCA) based on the compressed mixed linear model (MLM) were conducted using the R package of Genomic Association and Prediction Integrated Tool (GAPIT) (Lipka et al. 2013; R Core Team 2019). A PCA plot visualization was also performed in R software. LD decay analysis to identify candidate regions was performed using PopLDdecay version 3.27 (Zhang et al. 2019).

GWAS analysis

A GWAS analysis was performed to analyze the associations between genotype and phenotype using the GAPIT package (version 3.0) in R (Lipka et al. 2013; Liu et al. 2016). In addition, the MLM was adopted, and PCA and kinship were also used to validate population strati-fication with the GAPIT. To screen significant SNPs, the genotype was filtered by software PLINK (Purcell et al. 2007). Noneffective SNPs were removed, and a total of 97,469 effective and independent SNPs remained. A genome-wide threshold was calculated using the formula: “‒log10 (1 / number of effective SNPs)”. Finally, markers with an adjusted ‒log10 (P-value) ≥ 4 were regarded as the significant ones to identify the association loci and SNP markers located at locus peaks (Yang et al. 2014).

Haplotype analysis

For the haplotype analysis, SNP markers, except missing and heterozygote SNPs, were used to perform the analysis. The average score and variety count were determined from phenotype data for each variety, and haplotypes were identified that were significantly associated with the phenotype. The haplotype variation analysis was performed using PopART software (Leigh and Bryant 2015). The online tool Gene Structure Display Server 2.0 (Hu et al. 2015) was used to visualize the gene structure and SNP position.

Variation in survival rate under the submergence condition

The following six rice subgroups encompassing 137 rice accessions were used to assess submergence tolerance: Tropical japonica, temperate japonica, indica, aus, aromatic, and admixture. A wide range of submergence tolerance in different rice accessions was observed at the seedling stage under submergence condition. The survival rate ranged from 0 to 100% with an average of 40.58%. (Fig. 1A). The distribution of survival rate by subgroups was analyzed. The survival rate of Indica (43 accessions), Temperate japonica (62 accessions), Tropical japonica (19 accessions), Admixture (2 accessions), Aromatic (3 accessions), Aus (8 accessions) was 41.55, 42.69, 40.70, 36.67, 33.33, and 22.50, respectively. No significant difference was detected between subgroups by the Duncan test (Fig. 1B).

Figure 1. (A) Distribution of survival rate for the 137 rice accessions. (B) Box plot of total survival rate by subgroups. Adm: admixture, Aro: aromatic, Ind: indica, Tej: temperate japonica, Trj: tropical japonica.

Population structure and LD decay analysis

CV analysis was conducted using ADMIXTURE version 1.3.0, which indicated K = 6 was the optimal population group (Fig. 2A). The PCA analysis showed that the top two PCs each explained (61.86 and 25.12)%, which could explain most of the variation to select for visualiza-tion. Significant clusters belonging to the tropical japonica, temperate japonica, indica, and aus subgroups were observed in the PCA analysis (Fig. 2B). The population structure plot of 137 accessions was generated using Structure Plot V2.0, which was also divided into six groups that distinguished their subgroups. Temperate japonica was divided into Clusters 1 and 5, indica was divided into Clusters 3 and 4, admixture was divided into Clusters 2 and 3, and aromatic was divided into Cluster 2 and 6. Tropical japonica was dominant in Cluster 2, while aus was dominant in Cluster 6 (Fig. 2C). The estimate of genome- wide LD decay along physical distance was calculated using r² of allele pairs between two loci for the 137 rice accessions. The maximal r² value was 0.52, the threshold value was determined to be 0.26 half of the maximal r² value, and the LD decay distance was about 230 kb for a genomic candidate region (Supplementary Fig. S1).

Figure 2. Population structure analysis based on 137 rice accessions. (A) CV error of diverse groups (K). The dotted transverse line represents the lowest level. (B) Principal Component Analysis (PCA) (PC1 and PC2). red, brown, green, azure, blue, and purple represent the Adm, Aro, Aus, Ind, Tej, and Trj rice subgroups, respectively. (C) Plot for population structure analysis at K = 6.

GWAS for submergence tolerance

The GWAS analysis for submergence tolerance was conducted separately using the GAPIT package in R. The threshold was set as ‒log(p) ≥ 4 at a significance level of 0.01 after Bonferroni multiple test correction for signifi-cantly associated SNPs. In the Manhattan plots, nine lead SNPs were detected (Fig. 3). Considering the size of the LD block, the lead SNPs located inside the 460 kbp were regarded as being overlapped. Finally, nine lead SNPs were detected as QTLs for submergence tolerance as compared to the previously reported QTLs based on Gramene (http://archive.gramene.org, accessed on 1 November 2022) (Table 1).

Table 1 . QTL information detected in this study.

QTL	Lead SNP	Chr	‒LOG10(P)	Reported QTL		Reference of previously reported QTLs
				QTL ID	Related trait
qSUB3	26924139	3	4.42	qFlw3	Flag-leaf width	Mei et al. (2005)
qSUB4	22035756	4	4.11	orl1	Culm/leaf	Sato et al. (1998)
qSUB6-1	3321700	6	4.35	qbph6	Brown planthopper	Van mai et al. (2015)
qSUB6-2	26849774	6	4.05	qREP-6	Root elongation	Shimizu et al. (2008)
qSUB11-1	9322120	11	4.54	qDLR11	Alkaline stress	Qi et al. (2008)
qSUB11-2	18566175	11	4.45	yld11.1	Yield per plant	Moncada et al. (2001)
qSUB11-3	22317371	11	4.39	Xa21	Bacterial blight resistance	Ronald et al. (1992)
qSUB12-1	1392303	12	4.29	Pi-4(t)	Blast resistance	Inukai et al. (1994)
qSUB12-2	21636751	12	4.16	qSUB12.1	Submergence tolerance	Septiningsih et al. (2012)

Figure 3. The manhatton plot and QQ plot for GWAS analysis. (A) Manhattan plots for submergence tolerance in 137 rice varieties. X axis indicates physically mapped chromosomes. Y axis indicates significance as calculated by ‒log10 (p). (B) QQ plot for submergence tolerance.

Haplotype analysis of candidate genes

From the results of the haplotype analysis, two candidate genes showed significant differences among the groups of haplotypes. Finally, these two genes were detected as candidate genes for submergence tolerance. One of these candidate genes, Os03g0679300, contained four SNPs in the exon region (Fig. 4A). The haplotype analysis of the accessions showed that the four SNPs divided into four haplotypes, with the maximum phenotypic survival rate variation of 41.93 % between haplotype 2 and haplotype 4 (Hap 2 and Hap 4) (Fig. 4B). The rate of survival of plants was significant with Hap 1, Hap 2 and Hap 3 but was not significant with Hap 4. This was expected, as the major constituents of this haplotype were the aus and indica varieties (Fig. 4C). The Os11g0517800 gene contained 12 SNPs in the exon region (Fig. 5A). Twelve SNPs were divided into four haplotypes with maximum phenotypic variations of 72.22 % for survival rate between Hap 1 and Hap 3 (Fig. 5B). Hap 1 was the superior genotype in temperature japonica varieties (Fig. 5C).

Figure 4. Haplotype analysis of Os03g0679300. (A) Schematic representation of gene structure and SNPs positions in Os03g0679300. Yellow and blue blocks represent exon and untranslated region, respectively. Black vertical bars represent SNPs. (B) Results of haplotype analysis of Os03g0679300. (C) Haplotype variation analysis. Colors indicate rice subgroups as indicated in the legend. Circle size indicates the number of varieties in each Hap. Traverse lines represent the extent of variation between four haplotypes.

Figure 5. Haplotype analysis of Os11g0517800. (A) Schematic representation of gene structure and SNPs positions in Os11g0517800. Yellow block represents exon and black vertical bars represent SNPs. (B) Results of haplotype analysis of Os11g0517800. (C) Haplotype variation analysis. Colors indicate rice subgroups as indicated in the legend. Circle size indicates the number of varieties in each Hap. Traverse lines represent the extent of variation between four haplotypes.

Most of the low-lying areas and fragile areas with heavy rainfall in South and Southeast Asia are reported to be vulnerable to submergence, causing huge losses, and climate change is expected to worsen. In this study, we evaluated the phenotypic variations as submergence tolerance for the rice seedling stage.

We detected nine QTL regions for submergence tolerance based on GWAS analysis and haplotype analysis. Among the detected nine QTLs, three QTL regions (qSUB6-1, qSUB12-1, and qSUB12-2) were contained with previous reported genes associated with submergence tolerance. In the qSUB6-1 region, Os06g0164900 was reported to have S-Domain receptor-like kinase-66 (SDRLK-66) that responds to chilling and submergence. Also, in the qSUB12-1 region, Os12g0130500 was reported to have S-Domain receptor like protein-6 (SDRLP-6) that responds to submergence (Mondal et al. 2021; Naithani et al. 2021). One of the detected QTLs overlapped with the reported qSUB12.1 for submergence tolerance. This QTL was confirmed from the 466 F_2:3 population between IR72 and Madabaru (Septiningsih et al. 2012). In this QTL region, Os12g0541300 was associated with ethylene-induced aerenchyma forma-tion in roots under oxygen-deficient conditions (Yamauchi et al. 2017). Despite the haplotype analysis, we did not select a candidate gene in other QTL regions, except two QTLs that detected candidate genes.

One of the candidate genes, Os03g0679300, which is located in the QTL region of qSUB3 was reported to be similar to Ubiquitin-conjugating enzyme family protein; this protein is a key factor in ATP-dependent substrate proteolysis and is universally present in many organisms (Liu et al. 2020). Also, ubiquitin was reported to regulate various biological functions of the plant as diverse as abiotic stress response, plant growth, and development (Stone 2014; Shu and Yang 2017). Liu et al. 2021 reported that ubiquitin regulates the stability of the transcription factor to modulate the submergence response in the Arabidopsis. Os11g0517800 gene, located in the QTL of qSUB11-2, encodes the domain region of GCN2, which is regulated by the eIF-2-alpha kinase associated with the resistance of abiotic/biotic stress in rice. According to the latest study of Os11g0517800, it reported candidate gene of the bacterial blight resistance gene in rice (Park et al. 2022).

The haplotype analysis of the two candidate genes showed clear grouping by statistical analysis. Evaluations of these candidate genes might provide future strategies for developing submergence tolerant rice varieties.

In this study, the trait of submergence tolerance in rice was surveyed in a panel of 137 accessions, which identified nine QTL regions associated with submergence tolerance, and two candidate genes on chromosomes 3 and 11. Evaluations of the two candidate genes and QTLs reported in this study might provide strategies for future studies and breeding programs.

This work was supported for 2 years by Pusan National University.