Identiﬁcation of QTL for Tolerance to Flooding Stress at Seedling Stage of Soybean ( Glycine max L. Merr.)

: Flooding stress is a serious problem in soybean production, causing a remarkable yield reduction. The onset of rainy season during the early growth of soybean in Korea and some other parts of the world potentially subjects soybean plants to ﬂooding stress. The objective of this study was to map quantitative trait loci (QTL) for ﬂooding tolerance using a recombinant inbred line (RIL) population derived from a cross between ‘Danbaekkong’ (ﬂood-tolerant) and ‘NTS1116 (cid:48) (ﬂood-susceptible) cultivars grown in a plastic house for two years. The plants were ﬂood-stressed at the V1-V2 stage by ponding about 10 cm water from the soil surface. Leaf chlorophyll content and shoot dry weight were measured under control and ﬂooded conditions to map the QTL. The genetic map was constructed using 1689 polymorphic markers obtained from the 180K Axiom ® SoyaSNP markers used for genotyping the population. Ten QTL with 3.39–5.14 logarithm of odds scores and 8.1–30.7% phenotypic variations (PVE) were identiﬁed on seven chromosomes. One QTL on chromosomes 6 and 15 and two QTL on chromosome 7 were detected at least in two different environments causing up to 30.7% PVE, suggesting their potential applications in the breeding of ﬂood-tolerant soybeans. The results could be useful in further exploring the genetic basis of ﬂooding tolerance and developing tolerant cultivars


Introduction
Increments in the production of major crops like wheat, corn, rice, and soybean are important for global food security. However, the yields of many crops are negatively affected by various factors including climate change, which is supposed to aggravate several biotic and abiotic stresses, like flooding [1,2]. Various efforts, including the development of tolerant cultivars, have been made to increase the yield of soybean by minimizing the negative effects of those stresses [3][4][5]. Flooding stress accounts for a significant yield reduction in soybean [6,7] because soybean plants are sensitive to flooding stress at different stages, such as germination [8,9] and early vegetative as well as reproductive stages [10,11]. The coincidence of heavy rainfalls during the early growth stages makes the soybean plants liable to flooding stress in lowland, upland, as well as paddy fields under poor drainage conditions. In particular, the cultivation of upland crops in paddy fields is more challenging in some countries such as Korea and Japan [12,13]. The poorly drainable nature of paddy soils renders the soybean plants additionally prone to flooding stress. Hence, the development of flood-tolerant soybeans could be of great importance to increase soybean yields under flooding stress conditions. Flooding makes the soil environment hypoxic and the existence of prolonged flooding results in anoxic condition, retarding plant growth and development due to insufficient oxygen supply [14]. The retarded growth and development under the little or no oxygen environments is principally associated with the reduced ability of plants to metabolize nitrogen [15,16]. Although legume crops, like soybean, are capable of fixing atmospheric nitrogen into the soil, they require a sufficient amount of oxygen to accomplish such fixation, which is retarded under hypoxia and anoxia conditions [17]. The little or no oxygen conditions not only lower the nitrogen fixation but also reduce nitrogen uptake [18,19]. Additionally, a shortage of oxygen limits adenosine triphosphate synthesis, which upon conversion to adenosine diphosphate generates the energy required for plant growth and development [20]. The reduction in oxygen concentration under flooded conditions may result in elevated carbon dioxide concentration, creating a toxic soil environment [21,22]. The prevalence of a low-oxygen environment is far less detrimental than the carbon dioxide toxicity, which causes chlorosis, retards root growth, and may lead to the death of soybean plants [21].
The tolerance to flooding stress is a complex quantitative trait that is regulated by multiple genes, and is also greatly influenced by environmental conditions [23,24]. Identification of the genetic regions accounting for quantitative trait variations is useful to understand the genetic make-up of a trait and subsequently to improve the desired trait through marker-assisted selection (MAS). A few genetic mapping studies for flooding tolerance in soybean have been carried out [25][26][27][28]. Most of the studies are based on indirect measurements of phenotypic traits like yield variation or visual differences in the plant injury levels resulted from flooding stress [11,12,[25][26][27][28]. However, selection and accurate measurements of appropriate traits are crucial to precisely identifying quantitative trait loci (QTL) for stress tolerance [29,30]. Leaf chlorophyll content (CC) and shoot dry weight (DW) are among the traits that are highly influenced by flooding stress [25,31,32]. The CC, DW, and/or their index values have been considered for identifying the flood-tolerant QTL in soybean [33,34] as well as in other crops like common bean [35], wheat [36][37][38], and barley [39]. A similar study for identifying flooding tolerance-related QTL was conducted by Dhungana et al. [34] using a common parent "NTS1116" that was also used in the present study. It is worthwhile to contrast and compare the results of these studies because the results of QTL may be influenced by genetic background and/or environment [40,41]. Therefore, additional QTL information from comparable genetic backgrounds could be more useful in the flooding tolerance-related soybean breeding programs. The objective of this study was to identify QTL for flooding stress tolerance at the seedling stage of soybean using a recombinant inbred line (RIL) population.

Plant Materials
A RIL population was developed through the single seed descent method using 'Danbaekkong' (female parent) and 'NTS1116' (male parent) cultivars. The population consisting of 152 lines (F 6:7 and F 7:8 ) was used to map the QTL for flooding tolerance. The parents for developing the RIL population were selected based on a previous report [42], in which 192 soybean germplasms were screened for flooding tolerance and found that 'Danbaekkong' was a tolerant and 'NTS1116' was a susceptible cultivar. Furthermore, 'Danbaekkong' and the other flood-tolerant parent 'Paldalkong' that was used in a previous study [34] are among the elite soybean cultivar and parental line, respectively, in Korea [43]. Additionally, 'Danbaekkong' is resistant to bacterial leaf pustule [44] and 'Paldalkong' is resistant against soybean mosaic virus [45].

Phenotyping
The parents and RILs were grown in a plastic house under ambient temperature and light conditions in 2017 and 2018. In June of 2017, seeds were sown in stainless steel containers of 2 × 1 × 0.6 m dimensions at the 5-cm spacing between rows and plants.
Since the soil temperature in stainless steel and plastic containers (1.58 × 1.13 × 0.6 m) was not significantly different, both types of containers were used to grow plants in 2018. The possible undetectable variation in soil environments between the metal and plastic containers was further reduced by growing the same genotypes in the same type of containers under both flood-stressed and control conditions. The plant spacing was increased in 2018 to keep 10 cm from row to row and 7.5 cm from plant to plant, and seeds were sown in May. In both years, the containers were filled with the soil mixture composed of equal proportions of soil, sand, and compost. Plants were grown under normal conditions until the V1-V2, and then the stress-designated plants were exposed to flooding for 14 d by holding about 10 cm water from the soil surface. On the other hand, the control-designated plants were allowed to grow under normal conditions throughout the experiment period. In 2017, a single plant was considered as one replication; however, in 2018, an average value of the adjacently grown three plants of the same genotype was taken as one replication. In both years, plants were grown in a complete randomized design with three replications under both control and flood-stressed conditions.
The CC and DW data collected from the plants grown under the control and flooded conditions were considered to map the QTL. A chlorophyll meter (SPAD-502Plus, Minolta Camera Co., Osaka, Japan), which is widely considered for the determination of greenness of plant leaves, was used to measure the CC of soybean plants. The CC was measured on the middle leaflet of second and third trifoliate leaves at 12 and 13 days after flooding (DAF), respectively. In order to minimize the possible variation in CC due to soil plant analysis development (SPAD) reading, the two values obtained from two different leaves of each plant were averaged. For the DW measurement, plants were harvested at 14 DAF. Above-ground portions of the plants were kept into paper bags and oven-dried (60 • C) to constant weight. In addition to the CC and DW, their index values (CCI and DWI, respectively) were also calculated and included in the QTL analysis. The CCI and DWI were separately calculated as the ratio of their values (average of three replications) under flood-stressed to control conditions. Flooding tolerance index (FTI) was calculated as the mean value of CCI and DWI and was also considered for the QTL analysis.

DNA Extraction and Genotyping
For genomic DNA extraction, parents and RILs were grown in trays until the V2 stage and trifoliate leaves from 3-4 plants of each genotype were bulk harvested and stored at −80 • C. About~100 mg of the bulked leaves were ground in 2-mL tubes using a bead beater (TissueLyser II, Qiagen, Hilden, Germany). The powdered leaves were utilized to extract DNA using a kit (Exgene TM Plant SV Miniprep Kit, GeneAll, Seoul, Korea) following the manufacturer's instructions. The DNA was eluted with 50 µL Buffer AE supplied in the kit. The parents and 152 RILs were genotyped using 180K Axiom ® SoyaSNP array [46].

Data Analysis
Analysis of variation (ANOVA), construction of frequency distribution of the RIL population, and correlation analysis were done using SAS9.4 (SAS Institute Inc., Cary, NC, USA). Broad-sense heritability (H 2 ) was calculated using an equation described earlier [47] with some modifications.

Linkage Mapping and QTL Analysis
Only the markers that were polymorphic between two parents were selected from the 180K SNP. The polymorphic markers were further analyzed to find redundant markers. The redundant markers have identical segregation patterns and cluster at the same genetic position while constructing a linkage map and thus cannot provide an additional contribution to a genetic study [48]. Therefore, the redundant markers were separated using the Bin function in IciMapping [49] before a map construction. The criteria for removing the redundant markers using the Bin function were set as the significant segregation distortion of p < 0.001 and missing data with >10%. After removing the redundant markers, a linkage map was constructed using the Map function of IciMapping following the manufacturer's instructions. The parameters for the Map function were adjusted as grouping by 3.0 logarithm of odds (LOD) threshold, ordering by nnTwoOpt and rippling by the sum of adjacent recombination fractions.
The QTL for flooding tolerance were identified using composite interval mapping (CIM) in QTL Cartographer V2.5 (available at http://statgen.ncsu.edu/qtlcart/, accessed on 6 February 2020) following the manufacturer's instructions with adjusted parameters: Model 6, forward and backward regression, walk speed of 1.0 cM, and putative QTL with a window size of 10 cM. The number of control markers was 5, which was a default parameter. The LOD threshold for each trait was determined using a 1000 permutation test at p < 0.05. After the completion of the analysis, the QTL information was extracted by adjusting a minimum of 20 cM between QTL and 2-LOD support interval. The QTL were named following the method applied in an earlier report [34].

Screening of Candidate Genes
Candidate genes residing within the QTL regions were searched using SoyBase (www.soybase.org, accessed on 2 May 2021) and Phytozome (www.phytozome.net, accessed on 2 May 2021). Since chromosomes 3 and 7 harbored more than one QTL with at least one of >10% phenotypic variation (PVE), the candidate genes were searched in the QTL regions of these two chromosomes. The genes were searched based on their functions related to stress response, chlorophyll content, photosynthesis, and/or plant growth. The Glyma 1.1 gene version was used to collect the gene information.

Phenotypic Variation and Correlation between Traits
The CC of parents and RILs were higher in 2018 than that in 2017 under the control as well as flooded conditions; however, the DW were higher in 2017 (Table 1). Under the control condition, the CC of 'NTS1116' (33.38) was higher than that of 'Danbaekkong' (31.92)  The values of CCI, DWI, and FTI for 'Danbaekkong' were higher than those for 'NTS1116' in both years ( Table 2). The CCI, DWI, and FTI values of RILs were higher in 2017 (0.66, 0.79, and 0.72) than in 2018 (0.55, 0.47, and 0.51), respectively.
The broad-sense heritability (H 2 ) was calculated on a line-mean basis and found to be almost equal for both traits, CC (68.35%) and DW (68.42%) ( Table 1). The variations in CC and DW among RILs over two years showed continuous distribution with transgressive segregation (Figure 1). The ANOVA for CC and DW under flooding treatment indicated that both traits were significantly (p < 0.05) different among RILs, over two years (except for year-wise effect on CC), and their interactions (Supplementary Table S1).

Linkage Mapping
A total of 22,587 (12.52%) SNP markers were found polymorphic between the parents, among the 180,375 markers genotyped. The polymorphic markers were subjected to Bin function to remove the redundant markers (Supplementary Table S2), and thus 1689 SNPs remained for constructing linkage maps of 20 chromosomes (Supplementary Figure S1). The cumulative length of linkage maps was 4293.98 cM with an average of 2.54 cM between the markers. The longest (392.82 cM) and shortest (88.16 cM) linkage groups were formed by chromosomes 4 and 16, respectively.

QTL Analysis under Flooding and Index (CCI, DWI, and FTI)
Ten QTL for flooding tolerance were identified which were distributed on seven chromosomes (Supplementary Figure S2; Table 3). The QTL identified in different environments at the same, overlapping or adjacent markers were reported as the same QTL. They showed ranges of 3.39 to 5.14 LOD scores and 8.1 to 30.7% PVE.

Screening of Candidate Genes in the QTL Hotspots
The genomic site harboring a relatively large number of QTL within a region of a chromosome was considered as a QTL hot spot. Chromosomes 3 and 7 contained at least two major QTL (>10% PVE), and thus were chose to search the candidate genes. The functional annotation of the genes, resided in the QTL hot spots, is related to stress response, chlorophyll biosynthesis, and growth-regulating factor (Supplementary Table S3).

Discussion
Flooding stress negatively affects chlorophyll content and reduces biomass accumulation and seed yield in soybean [32,50]. Development of flooding tolerant soybean cultivars could be of great significance to minimize yield loss, for which the collection of information on genetic basis of flooding tolerance through QTL analysis would be highly useful.
The H 2 for CC and DW might have been influenced by the significant interaction of RIL × year, indicating substantial effects of the environment on CC and DW. The difference in air temperatures between two years might have played a role for the significant RIL × year interaction. Heritability can change over time because of the change in genetic variance due to environmental factors that may affect the interaction between genes and the environment [51]. The mean monthly temperature in 2017 (32.2 • C) was higher than that in 2018 (28.7 • C) during the plant growing period. The higher CC in 2018 was possibly owing to the lower temperatures [52]. The inconsistent correlations observed between CC and DW under control and flooded conditions over two years was also reflected in higher DW value of RILs under the control but lower value under flooding stress in 2018 than in 2017. Similar results were observed in our previous study where RILs of mapping population with different genetic background were phenotyped for flooding tolerance in an early growth stage [34].
Distribution of the RIL population for CC and DW values showed a continuous variation in the traits with transgressive segregation, which could be utilized to select RILs having extreme phenotypes. Based on the mean FTI value of two years, 41 and 10 RILs showed positive and negative transgressive segregation, respectively (Supplementary Table  S4). Furthermore, the transgressive segregation obtained in the present study may provide novel soybean lines that could be exploited in flooding tolerance breeding programs [34,53]. Differences in CCI, DWI, and FTI between two parents and wide ranges of the values in the RIL population in both years depicted an appropriateness of the population in mapping QTL for flooding tolerance. An adoption of quantitative phenotyping method instead of visual observations to measure the flooding tolerance could increase the reliability of the QTL results.
Five QTL identified in this study overlapped the chromosomal locations or mapped in the vicinity of candidate genes and/or QTL associated with tolerance to abiotic stresses, including flooding. A major QTL qSFT_7-3 (up to 30.7% PVE) co-localized with a previously reported QTL for flooding tolerance in soybean [54]. These major QTL cover several candidate genes (DEAD-box ATP-dependent RNA helicase 31, WRKY transcription factor 55, Dof domain, zinc finger protein, and B-box zinc finger protein) associated with abiotic stress in different plants, including soybean [55][56][57][58][59][60]. Reports show that some genes respond to multiple abiotic stresses, for example VlWRKY3 gene promotes salt and drought stress tolerance during the germination, seedling, and the mature plant stages of transgenic Arabidopsis thaliana [61] core abiotic stress-responsive genes enhance drought, waterlogging, and osmotic stresses tolerance in sesame [62], and the DGK gene family responds under polyethylene glycol, salt, alkali, and salt/alkali treatments in soybean [63]. Candidate genes Glyma03g34900, Glyma03g35010, and Glyma03g35180 related to stress and/or growthregulating factor are found in the QTL qSFT_3-64. The QTL on chromosome 7 harbor several candidate genes associated with root and shoot growth factor, chlorophyll fluorescence, stress response, light-harvesting, and chlorophyll a/b binding protein (Supplementary  Table S3), indicating their roles in CC and DW of soybean plants. The QTL qSFT_7-3 was also found to cover the chromosomal location of SNP markers (causing up to 5.7% PVE) associated with flooding tolerance at an early reproductive stage in soybean [54]. A QTL causing up to 10.9% PVE for flooding yield index in soybean [11] was found within the QTL qSFT_13-53. Few stress tolerance-related candidate genes like Glyma13g27120, Glyma13g27410, and Glyma13g29760 also reside in the QTL regions of chromosome 13. A candidate gene Glyma15g10750 encoding for methionine sulfoxide reductase was found within QTL qSFT_15-67, implying its role in flooding tolerance. The methionine sulfoxide reductase is responsive to abiotic stresses in soybean [64][65][66].
QTL qSFT_7-3 identified in the QTL hot spots, overlapped the physical positions of the QTL detected in previous studies. In addition, two SNPs (AX-90363591 and AX-90445362) on chromosome 13 that were linked with flooding tolerance QTL [34] were also found to be linked with the QTL in this study (Supplementary Figure S3). Furthermore, an RNA-sequencing report [67] revealed that a few abiotic stress-related genes xyloglucan endotransglucosylase/hydrolase (LOC100799976), polygalacturonase (GLYMA_07G066900), calmodulin-binding protein (GLYMA_07G008400), soybean NAC gene (GLYMA_07G050600), DEAD-box RNA helicase gene (GLYMA_07G056600), and dehydration-responsive element-binding (GLYMA_07G017300), which were upregulated in the tolerant parent 'Danbaekkong' and downregulated in the susceptible parent 'NTS1116', reside in the QTL hot spots on chromosome 7. This RNA-sequencing report also supports the reliability of the QTL identified in the present study. The co-localization of flooding tolerance-related QTL identified in the present and previous studies, conducted at the same and/or different stages of soybean growth using diverse genetic background, indicated that the QTL could be of practical application and be utilized in improving the flooding tolerance of soybean cultivars. Although some of the QTL identified in this study physically located in or near the QTL or candidate genes for flooding tolerance, many of them were not consistently identified over different environments possibly because of the involvement of several biochemical mechanisms in abiotic stress tolerance [68] and influence of growing environment [69]. Moreover, several genes of different metabolic pathways like SnRK1A [70], N-end-rule [71,72], and glycolytic [73] are reported to interact for flooding tolerance.

Conclusions
A total of 10 QTL for flooding tolerance were identified on seven chromosomes. Some QTL were found to overlap while others to reside in the vicinity of some stress-related QTL and/or candidate genes reported earlier. At least one QTL on chromosomes 6 and 7 were identified in more than one environment. Since the abiotic stress tolerance in plants is a complicated and multi-dimensional phenomenon, analysis of the QTL for flooding tolerance by using the derived phenotypic data (such as CCI, DWI, and FTI) might be as equally useful as the data obtained from direct measurements (for example, CC and DW). This study is applicable in understanding the genetic basis of flooding tolerance and in utilizing the QTL in breeding programs. Research is under way to analyze the expression of candidate genes residing in the QTL hot spots, using flooding stress tolerant and susceptible RILs grown under control and flooded conditions. The QTL could be transferred into soybean genotypes of interest to enhance their flooding tolerance by adopting MAS technology. The application of this study may even increase in the context of climate change and potential consequences of increased flood events.  Table S2: Redundant SNP markers removed by the Bin function of ICI mapping; Table S3: Candidate genes residing in the QTL hot spots on chromosomes 3 and 7;