DNA was isolated from the leaves of two parents and 186 F₁₀ RILs at the tillering stage by the CTAB method. DNA quality was detected by agarose gel electrophoresis and nanodrop nd1000 spectrophotometer. In this study, the genome re-sequencing depth of parents TD70 and Kasalath is ~20×. The average sequencing depth of RILs is 10×. Sequencing was carried out with Illumina sequencing platform NovaSeq 6000, and the sequencing mode is paid end 150 bp. After low-quality base filtration of sequenced reads, they were compared to the reference genome using BWA software (Li and Durbin, 2009) (http://rapdb.dna.affrc.go.jp/ , IRGSP-1.0 version ). Then, SAMtools (Li et al., 2009) and PICARD (http://picard.sourceforge.net ) were used to remove redundancy. GATK (McKenna et al, 2010) was performed for calling SNP and genotype identification.

The resequencing of 186 RIL lines from TD70/Kasalath was performed for SNP identification, and high-quality, biallelic homozygous loci were obtained. To avoid false positive SNP genotyping in the population, Huang et al. (2009) used the sliding window method and made appropriate modifications to evaluate a group of continuous SNPs and perform genotyping. The genotype of RILs is split into a series of recombinant bins according to recombination breakpoints. The sliding-window SNP number of these bins is 15 and adjacent Bins within 200 kb are merged. The bin genotype of the transition region between two different genotype blocks is set as the missing data. The abnormal separation mark showed that the partial separation (P<0.01) was discarded.

QTL IciMapping v4.1 software (Meng et al., 2015) was used to construct the genetic map using the Kosambi map method and analyze the four traits related to coleoptiles on the 4-6th days after anaerobic stress treatment of the recombinant inbred lines. The inclusive composite interval mapping (ICIM-ADD) method was used for QTL identification. The logarithm of odds (LOD) score threshold was set to 2.5. Loci with LOD ≥ 2.5 are deemed candidate QTLs associated with the trait. QTL naming refers to McCuch’s method (McCouch, 2008).

In order to obtain expression profiles and differentially expressed genes related to anaerobic germination, the seeds of parents TD70 and Kasalath were taken for RNA Seq experiments with/without anaerobic treatment, respectively. The seeds were rinsed with double-distilled water, and then 5-6 seeds were taken from each sample, and three biological replicates were set, and thus a total of 12 samples were taken. These samples were immediately put into liquid nitrogen for quick freezing and stored at - 80 °C for RNA isolation. The purity, concentration and integrity of the total RNA isolated were evaluated, and RNA seq library construction and sequencing were referred to the previous process (Liang et al., 2022).

The Fastp program (Chen et al., 2018) was engaged to remove the adapter, low-quality sequencing tail, and then the Hiast2 program (Pertea et al., 2016) was used for mapping reads to the reference genome (http://rapdb.dna.affrc.go.jp/ , IRGSP-1.0 version). SAMtools (Li et al., 2009) were used for converting format. The expression abundance of transcripts was calculated based on fragments per kilobase of transcript per million fragments mapped (FPKM). The transcripts with a maximum FPKM value of less than 1 in multiple samples were eliminated as no expression or false positive transcripts (Trapnell et al., 2012). Differential expression genes were screened by the DESeq2 software package (Anders and Huber, 2010), with screening criteria of Log2 (fold change) ≥ 1 and q-value<0.05.

The faster the rice coleoptiles elongate, the sooner the seedlings escape from an anoxic environment and increase the survival chances of rice. Therefore, rice coleoptile is a classic tissue for studying AG tolerance (Narsai et al., 2015). In this study, we investigated the dynamic changes of coleoptiles of two parents under anaerobic conditions; Significant variation in coleoptile traits was observed between TD70 and Kasalath, which were the parents of both ( Figure 1 ).

After one day of treatment, Kasalath seeds began to germinate, and obvious coleoptiles bulged, while TD70 showed no significant change. On the second day, Kasalath showed obviously elongated coleoptiles, while TD70 seeds had just begun to germinate, and prominent coleoptiles emerged, which showed that Kasalath seeds germinated faster than those of TD70. The development of Kasalath coleoptile was occurred mainly from the second to the fifth day of treatment, while the growth of TD70 coleoptile mainly occurred on the third day. On the second day of treatment, the growth rate of coleoptile length (CL) of Kasalath was significantly higher than that of TD70, and the difference was very significant; However, the CL of TD70 increased rapidly from the third day to the fourth day of the treatment, and there was no significant difference between them on the fourth day. Then, on the fifth and sixth days, the CL length of kasalath was significantly longer than that of TD70, and the difference was significant. It could be seen that, on average, the CL elongation speed of Kasalath was faster ( Figures 1A, B ). From the third day of treatment, the CSA and CV of Kasalath were significantly higher than those of TD70. However, from the third day, the CD of TD70 was significantly larger than that of Kasalath, and the difference was significant or extremely significant until the sixth day.

Thus, the phenomenon stated above indicated that the coleoptile development of Kasalath mainly occurred on the first and second days of treatment. In contrast, TD70 mainly occurred on the third day ( Figure 1 ). Thus, the four agronomy traits with significant or extremely significant differences in AG anaerobic treatment between the parents TD70 and Kasalath indicate that there were great genetic differences between them, on which QTLs mapping was conductive based.

The four traits of the coleoptile phenotype showed different segregation in the RIL population ( Figure 2 ). The coefficient of variation (CV) on the fourth day of treatment ranged from 18.92% to 66.19%; The coefficient of variation at six days was 17.09% - 52.21%. Moreover, these traits show transgressive segregation in variated degrees. All characters are continuous unimodal distribution, and the absolute skewness value is close to 0 in each character, indicating that all these characters conform to the normal distribution and are controlled by multiple micro genes ( Table 1 ).

This study obtained 1,304 Gb high-quality sequencing data from illumina paired-end sequencing. The average sequencing depth is 18.8×. Among them, ~97.25% of the reads could be uniquely mapped to the reference genome. These SNP sites of RIL individuals were identified, and SNPs with depth >5 were retained. After filtering exceptions, 1,344,770 high-quality SNPs were finally used for subsequent analysis of recombination events and construct binmap, and the map contains 12,328 bin markers. All bin markers are evenly distributed on 12 chromosomes. The number of markers on chromosome 1 was the largest, 148,582; chromosome 9 has the least markers, 87,301( Table S2 ).

Using the inclusive composite interval mapping (ICIM) method, 50 QTL loci associated with CL, CD, CSA and CV were identified in two stages of anaerobic treatment. These QTLs are distributed on all chromosomes except chromosome 9. There were as many as 10 loci distributed on chromosome 4, followed by 7 loci detected on chromosome 8 and 12, respectively; only one epistasis was detected on chromosome 2. Among the four traits related to coleoptile, the number of QTLs detected by coleoptile volume was the largest (17), followed by coleoptile diameter (16), and there were 5 QTLs associated with coleoptile length. These QTLs can explain phenotypic contribution rates ranging from 0.34% to 11.17% and LOD values ranging from 2.52 to 11.57. Some QTL loci in this study showed a positive additive effect and some negative additive effect, suggesting that both parents contributed favorable alleles.

There are 5 QTLs associated with CL. The contribution rate of a single QTL ranged from 4.85% - 10.99%. There were 12 loci associated with CSA, and the phenotypic contribution rate was between 1.40% to 9.39%. Sixteen QTLs were associated with CD, which could explain the phenotypic variation of 2.05% - 11.17%. Seventeen QTLs were associated with CVs, with phenotypic contribution rates ranging from 0.34% to 9.37%. QTLs showing physical-position overlap were considered as mere one QTL site. Finally, 42 QTLs were obtained in two stages of anaerobic germination. It was remarkable that five QTLs were detected in two germination stages: qCD5, qCD8, qCV1-1, qCD10-2 and qCV8-1. Three loci were detected simultaneously by two different traits, of which qCSA6 and qCV6-5 overlapped as just one locus, and as well as qCV4-2 and qCL4-2, qCSA8 and qCV8-1, respectively; and they affected CL, CSA and CV ( Table 2 ; Figure 3 ).

We found 6 out of the 42 loci have been previously reported ( Table 2 ). The location of qCV1-2 overlaps with the physical interval of which qAGP1 (Liu et al., 2021). QCL4-2/qCV4-2 and qCV8-2 are in the same physical regions as qSFW-4 and qRD-8, respectively (Yang et al., 2021). The physical area of qCSA4-2 is consistent with that of qCL-4-1; The physical location of qCD12-3 coincides with that of qAG12/qCSA-12-2 (Hsu and Tung, 2015; Yang et al., 2019). qCSA7 locates on chromosome 7, and its physical location in a region from where many QTLs related to anaerobic germination were found (Ling et al., 2004; Angaji et al., 2010; Septiningsih et al., 2013; Hsu and Tung, 2015). The results of the above report indicate the accuracy of our mapping result.

In the present study, three kinds of QTLs were defined as stable QTLs, including overlapping QTLs associated with multiple traits; the loci detected in both periods and the loci colocalized with reported QTL loci. There were thus 12 stable loci discovered in this study, and they were associated with 20 QTLs ( Table 3 ). To further clarify the effects of these QTLs, we summarized the phenotypic differences between two alleles at each locus in the RIL population.

Firstly, the RIL population was divided into TD70 type and Kasalath type according to the genotype of each marker site, and then the corresponding traits were compared between the individuals of the two types. The result showed that these 8 of 12 loci significantly affected coleoptile traits. The average value of the traits in the individuals with excellent alleles was higher than in those without excellent alleles. For instance, qCL4-2/qCV4-2, the maximum absolute difference in coleoptile length between the two alleles was 0.36 cm, and the volume difference was 2.44 mm³ which was highly significant. The maximum difference in coleoptile volume of qCV1-2 between the two alleles was 1.91 mm³, and the maximum absolute difference of coleoptile length was 0.31 cm, which was significant. This shows that these 8 loci were reliable and responsible for phenotypic variation. The other four loci were associated with CV and CSA. The average phenotypes of individuals with excellent alleles were higher than those of individuals without ( Supplementary Table S3 ).

In order to further confirm the result shown above, the effects of these excellent alleles were evaluated by polymerization of different loci. If the interaction between these 12 stable loci were not considered, the more excellent alleles aggregated in one individual would result in a better phenotypic value of these traits. In this study, most of the lines in the RIL population contain 3-7 favorable alleles, and their corresponding phenotypic value increased along with the aggregation of favorable alleles ( Figure 4 ). In this study, 12 lines containing 7 favorable alleles were obtained. Under hypoxia conditions, the length of coleoptiles grew rapidly, and their surface area and volume were also larger. Hence their AG tolerance was strong ( Figure 4 ; Supplementary Table S4 ). These results indicate that aggregation of favorable alleles could improve AG tolerance.

For identifying candidate genes associated with AG traits, we further investigated transcriptome response during the process of tolerance to anaerobic germination. As a result, 2,405 DEGs from kasalath under AG treatment were identified and consider as Kas-specific DEGs response to AG. Similarly, 2,620 TD70-specific AG-responsive genes were also identified ( Figure 5D ). GO and KEGG enrichment were further employed for dissecting the potential molecular mechanism under AG treatment. From GO/KEGG enrichment results, amount of DEGs were found to be invovled in many biological processes such as “response to oxidative stress, GO:0006979”, “hormone-medidated signaling pathway, GO:0009755”, “detoxification, GO:0098754”, “Amino sugar and nucleotide sugar metabolism” ( Figure 5 ). Combined with 12 stable QTLs, 103 expressed genes were located in QTL regions, of which 31 genes such as OsUsp1, OsMADS4, PR5, OsSPARK4, OsSUAR56, TOND1 show significant among different groups ( Supplementary Table S6 ).

Rice varieties with strong tolerance to anaerobic germination (AG) are suitable for direct sowing. So far, there are two methods to identify rice AG tolerance. One is according to the survival rate of 10-20 cm under submerged treatment for 21 days, which has been used to identify several QTLs (Angaji et al., 2010). However, this method is time-consuming and laborious. Another is to measure the coleoptile length. The QTL associated with coleoptile length coincided with that detected by using survival rate (Hsu and Tung, 2015). In order to study rice AG tolerance more comprehensively, four traits (CL, CD, CSA, CV) relevant to coleoptile were analyzed in this study.

Firstly, we identified the coleoptile phenotype of two parental materials under anaerobic treatment. The result showed that Kasalath germinated earlier than TD70, and the coleoptile grew rapidly. On the sixth day of treatment, the length, surface area and volume of Kasalath coleoptile were significantly higher than TD70 (but coleoptile diameter), which indicated that Kasalath had stronger AG tolerance than TD70 ( Figure 1 ). Further analysis of the four traits of coleoptile in RIL populations under AG conditions showed that they segregated into various degrees, which showed the characteristics of normal distribution. Therefore, the result indicated it was possible to use Kasalath and TD70 construction population to locate QTL of AG tolerance.

Recently, the next-generation sequencing technology (NGS) has been used for population genotyping, which is improving the accuracy of QTL mapping (Huang et al., 2009; Wang et al., 2011). In this study, we sequenced 186 RILs and detected 1,344,770 high-quality population SNPs with uniform distribution in the whole genome. The recombination breakpoint was determined by checking the position of genotype change. This method transformed the original SNP markers into effective recombination bins, which could be considered an effective genetic marker (Wang et al., 2011). With this method, a high-density genetic map containing 12,328 bin markers was constructed in this study, and the mapping of QTLs was conducted using this map.

In this study, we measured four coleoptile traits, CL, CD, CSA and CV, on the fourth and sixth days of anaerobic treatment. The results showed that the phenotypic variations of CL, CD, CSA and CV were normal distribution in RILs, indicating that they were consistent with the characteristics of quantitative genetics controlled by multiple genes ( Figure 2 ). In this population, 50 QTLs associated with coleoptile traits under the anaerobic were identified, and 12 stable loci, including 20 QTLs, were obtained through further analysis.

Hypoxia under anaerobic conditions changes plant metabolism, thus affecting plant growth. The change in the glycolysis pathway is one of the typical metabolomic responses to hypoxia. Upregulation of amylase, trehalose-6-phosphate phosphatase, hexokinase and other enzyme activities under hypoxia conditions provided sufficient evidence (Magneschi and Perata, 2009; Kretzschmar et al., 2015; Loreti et al., 2018). When the oxygen concentration is low, carbohydrate metabolism is subsequently inhibited, especially the step of starch decomposition into monosaccharides for glycolysis (Miro and Ismail, 2013). These studies identified complex mechanisms related to coleoptile growth, including carbohydrate metabolism, fermentation, hormone induction, cell division and expansion. Many transcriptome analyses are conducted to explore the molecular mechanism in regulating rice coleoptile growth under hypoxia and anaerobic conditions. The results showed that starch degradation and fermentation-related genes were up-regulated in anoxic rice seedlings (Narsai et al., 2015; Hsu and Tung, 2017).

Research shows that the elongation of rice coleoptiles is relevant to the hexokinase gene (HXK6, LOC_Os01g53930) under anaerobic conditions (Hsu and Tung, 2015). LOC_ Os06g03520, a candidate gene encoding DUF domain protein, is highly expressed in the coleoptile under flooded conditions (Zhang et al., 2017). OsTPP7 gene, encodes a trehalose-6-phosphate phosphatase. It indicates a low sugar utilization rate by increasing the turnover of trehalose-6-phosphate, which can also improve the pool capacity in heterotrophic tissue. Furthermore, it can enhance starch utilization, drive germinating embryo growth, prolong coleoptiles, and enhance anaerobic germination tolerance (Kretzschmar et al., 2015). AG tolerance can be achieved by combining the number of favorable alleles in different varieties.

Recently, OsGF14h, a key gene determining its strong germination and emergence under hypoxia, was cloned from japonica weedy rice via genome-wide associating study and QTL mapping. This gene encodes a “14-3-3 protein”, which is the upstream switch of abscisic acid (ABA) signal transduction. By interacting with two transcription factors, OsHOX3 and OsVP1, the gene inhibits the activity of the ABA receptor OsPYL5 in a hypoxic environment, reduces the sensitivity of ABA, and also activates the biosynthesis of gibberellin (GA), ultimately causing strong germination and emergence of weedy rice seeds. However, in modern genetic improvement, the lost of 4 bases of OsGF14h in cultivated japonica rice caused the the weakening of its function, so germination and emergence were limited under hypoxia conditions (Sun et al., 2022).