Genetic analysis of rice seedling traits related to machine transplanting under different seeding densities

Background Due to the diversity of rice varieties and cropping systems in China, the limitation of seeding density and seedling quality makes it hard to improve machine-transplanted efficiency. Previous studies have shown that indica and japonica varieties varied in machine transplanting efficiency and optimal seeding density. In this study, a RIL population derived from ‘9311’ and ‘Nipponbare’ were performed to explore the seedling traits variations and the genetic mechanism under three seeding densities. Results The parents and RIL population exhibited similar trends as the seeding density increased, including seedling height and first leaf sheath length increases, shoot dry weight and root dry weight decreases. Among the 37 QTLs for six traits detected under the three seeding densities, 12 QTLs were detected in both three seeding densities. Five QTL hotspots identified clustered within genomic regions on chromosomes 1, 2, 4, 6 and 11. Specific QTLs such as qRDW1.1 and qFLSL5.1 were detected under low and high seeding densities, respectively. Detailed analysis the QTL regions identified under specific seeding densities revealed several candidate genes involved in phytohormones signals and abiotic stress responses. Whole-genome additive effects showed that ‘9311’ contributed more loci enhancing trait performances than ‘Nipponbare’, indicating ‘9311’ was more sensitive to the seeding density than ‘Nipponbare’. The prevalence of negative epistasis effects indicated that the complementary two-locus homozygotes may not have marginal advantages over the means of the two parental genotypes. Conclusions Our results revealed the differences between indica rice and japonica rice seedling traits in response to seeding density. Several QTL hotspots involved in different traits and specific QTLs (such as qRDW1.1 and qFLSL5.1) in diverse seeding densities had been detected. Genome-wide additive and two-locus epistasis suggested a dynamic of the genetic control underlying different seeding densities. It was concluded that novel QTLs, additive and epistasis effects under specific seeding density would provide adequate information for rice seedling improvement during machine transplanting. Supplementary Information The online version contains supplementary material available at 10.1186/s12863-020-00952-1.


Background
Along with the social and economic development, rural labor transfer and ageing, traditional rice planting by hand technology has been unable meet the demands for rice production in China. Mechanization of rice cultivation is of great significance to improve rice production capacity and ensure national food security [1,2]. Until 2012, the comprehensive mechanization level of rice production in China was 68.82%, of which the level of machine tillage and harvest were 93.29 and 73.35%, respectively. However, the level of machine transplanting only accounted for 31.67%, which was the bottleneck of rice production mechanization [3]. Although rice mechanical transplanting technology has been explored since the 1950s in China, there are still many problems on cropping system and varieties application. The rice mechanical transplanting technology including mechanical equipment and seedling raising technology widely applied in northern China was derived from Japan based on conventional japonica rice, which was restricted applied in indica varieties and hybrid rice with double season and multiple cropping system in southern China [4][5][6][7][8][9]. In particular, due to the tight growing season, the limitation of seeding density and seedling quality makes it more difficult to improve the efficiency for mechanized transplanting [10][11][12][13]. Although it has been used in hybrid rice, which accounts for about 60% of the rice producing area, the traditional rice mechanical transplanting technology cannot fully exploit the high yield advantage of hybrid rice, mainly due to the poor seedling cultivation, such as high seeding density, poor quality of seedlings, high rate of seedlings injury and large amount of seedlings per hill [14][15][16][17]. Therefore, due to the diversity of growth characteristics of different types of cultivars, it is necessary to make the seedlings with uniform size, consistent growth, and seedling characteristics suitable for appropriate ecological zones, planting systems and planting methods under machine transplanting.
The most important factor limiting the efficiency of machine transplanting was seedling cultivation, and the optimum rice seedlings were restricted by some factors including cultivation methods, nursery substrates, variety types and seeding density [8,15,[18][19][20]. In general, it takes less time to recover from transplanting shock for young seedling with three to four leaves, seedling height ranging from 12 to 17 cm. Meanwhile, strong root crosslinked to stabilize the seedbed on the seedling-nursery tray not only reduced seedling damage during machine transplanting, but also conducive to the occurrence of tillers and the formation of yield [21]. Previous studies have shown that different rice varieties varied in machine transplanting efficiency and optimal seeding density [7,20,22]. Recently, the bowl-shaped blanket-like seedling transplanting technology combining the special machine transplanting seedling substrate, effectively solved the problems existing in the traditional blanket and potted seedling cultivation, including poor quality and quantitative positioning, high rate of seedling leakage and injury and reducing other adverse effects and the incidence of blight [1,23,24]. Seeding density was one of the most important constrains to produce good quality of rice seedlings in the seedling-nursery tray and were related to the characteristics of rice varieties [25]. Previous reports have suggested that seedlings with high density always have evenly emergence than those with the low density, and the consolidation force of the root system could be well formed into a blanket with the increase of seeding density. However, the quality of seedlings became worse with the increase of seeding density, and the seedling quality indexes such as dry matter of seedlings and root vigor showed a decreasing trend with the increase of seeding density [7,16,18,20,26]. Therefore, it is important to improve the efficiency of mechanical transplanting by selecting suitable seeding density for the special variety and raising method.
Many efforts have been made to improve seedling cultivation by providing plant growth regulators, which might cause several undesirable consequences such as seedling dwarfing, growth retardation and even reduction in yield if not used properly [27]. Therefore, it is urgent to investigate the genetic basis of seedling traits related to machine transplanting and improve the seedling quality through marker-assisted breeding. The essential problem of rice seedling response to seeding density was shade avoidance response, which referred to a set of architectural responses including accelerated growth of hypocotyl, internode, and petiole, decreased leaf surface area and chlorophyll and changed of leaf angle [28]. Competing for light and nutrients between adjacent plants at different seeding densities usually associated with balancing the investments in roots and leaves for rice seedlings. The molecular mechanisms of the hypocotyl, internode, or petiole elongation in shade avoidance response had been mainly characterized in dicots, depending on the cascade reaction of the light signal system, plant hormone signaling pathways, and growth regulation [29][30][31]. Previous studies showed that shading treatment through red to far-red light ratio (R/ FR) could induce rice seedling stem elongation, which was contributed by phytohormones signals [32]. Beyond that, planting density may also affect rice root system growth through root-root recognition mediated by root exudates [33]. Identification of quantitative trait loci (QTLs) and candidate genes associated with early vigor such as mesocotyl elongation, shoot length, stem length, fresh weight, dry weight, and root vigor had been undertaken using different types of mapping population in rice under direct-seeding system [34][35][36][37][38][39][40][41]. Although young seedling traits associated with dry weight and root vigor had been extensively researched under controlled laboratory conditions [25,42,43], there have been very few reports on combining genetic research for early seedling traits related to mechanical transplanting directly on seedling nursery tray as well as controlled condition such as seeding density or temperature [44,45]. The combined analysis would help in detecting QTLs that across different seeding densities as well as the special QTLs for certain traits that were important for mechanical transplanting seedling cultivation system. Therefore, the objectives of our research were applying a RIL population derived from indica and japonica cross: (1) to explore the trait performance across different seeding densities; (2) to compare the differences of QTLs for seedling traits under different seeding densities; (3) to analyze the whole-genome additive effects and epistatic effects under different seeding densities; (4) to detect novel loci controlling seedling traits under different seeding densities. These results would provide useful information for machine-transplanted rice seedling cultivation improvement in southern China rice-growing districts.

Seedling trait performance of RILs and their parents
The difference in seedling traits between the two parental cultivars ('9311' and 'NIP') was largely affected by seedling density (Additional file 1: Table S2 and Table 1). The two parents differed significantly at P ≤ 0.05 for all the six traits under three seeding densities. However, the degree of difference decreased when SLL and RDW were measured as the increasing of seeding density, while the degree of difference between the parents for SH, FLSL, FLL and SDW increased with increasing density. These results indicated the effect of seeding density on seedling traits was not consistent, thus suggesting significant differences of genetic mechanisms for seedling traits between inter-subspecies. In the RIL population, all the seedling traits showed continuous variation and obvious transgressive segregation, following approximate normal distributions. With the increase of seeding density, the average values of RIL population for SH, FLSL, FLL and SLL gradually increased, while the RDW and SDW gradually decreased (Fig. 1). The broad-sense heritabilities ranged from 71.37 to 90.91%, indicating the complexity of the genotypic response to seeding density. The heritability of FLL and SLL were moderate while it was low for RDW and SDW. There was no significant difference in heritability between different densities. The G × E interactions were highly significant (P ≤ 0.01) among the three seeding densities, suggesting the effect of seeding density on seedling traits should not be ignored ( Table 1). The correlations among different seedling traits showing similar trends between different seeding densities (Additional file 1: Table S3). FLL and SLL showed the highest correlations in both experiments. RDW was negatively correlated with other traits, but it was positively correlated with SDW. Taken together, these results indicated that seedling quality of the parents and RIL population on plastics nursery tray generally declined with the increase of seeding density, and the influences of seeding density on indica and japonica varied on traits. However, little was known about the genetic loci controlling seedling traits that were affected by seeding density, and the distributions of loci that differ in response to seeding density between indica and japonica rice genomes.
QTL analysis of seedling traits under three seeding densities QTL analysis based on SNP bin map under three seeding densities in three experiments identified 37 QTLs. Four QTLs for SH, 9 QTLs for FLSL, 6 QTLs for FLL, 5 QTLs for SLL, 5 QTLs for RDW and 8 QTLs for SDW were identified above the LOD thresholds under three seeding densities (Additional file 1: Table S4). We further analyzed in detail the QTLs that detected in at least two experiments ( Table 2).
QTL hotspot is a region where multiple traits were colocated. There were five QTL hotspots (bin311, bin411, bin1135, bin1565 and bin2518) identified clustered within genomic regions on chromosomes 1, 2, 4, 6 and 11 ( Table 2 and Fig. 2). However, it was not clear whether these five QTL hotspots were single locus with pleiotropic effects on the multiple traits, or a group of tightly linked loci. Further analyses revealed some noteworthy features about these regions. Three of these hotspots, bin311, bin1135 and bin2518 showed almost consistent effects across the three seeding densities, while bin411 mainly detected at MD and HD for FLL and SLL, with minor effects at LD. However, bin411 functioned at LD for SDW. The QTL hotspot at bin1565 on chromosome 6 were simultaneously detected for FLL and SLL at all three seeding densities, however, this hotspot was only detected for FLSL at LD and MD. These results suggest that bin311, bin1135 and bin2518 were not sensitive to seeding density, while bin411 and bin1565 appeared to be QTL-specific for seedling density. Another interesting finding was that the effects of QTL hotspots for the multiple seedling traits showed diversity. For bin311, the allele from 'NIP' increased trait values, with the additive effects in the same direction for SH and FLSL. For the other three QTL hotspots (bin1135, bin1565 and bin2518), the allele from '9311' increased trait values, and consistent among those multiple traits. However, the additive effect of bin411 showed the opposite direction, with positive effects for Fig. 1 The comparisons of seedling traits performances under three seeding densities. (a-f) Beanplot and a density plot with symmetrical arrangement of each traits, in which each blue beanline represented each observation value, red line represented the mean value. LD, MD, and HD indicated low, medium, and high seeding density, respectively. Ex1, Ex2 and Ex3 indicated the three experiments, respectively SDW and negative effects for FLL and SLL. All the five QTL hotspots produced considerable individual effects on the seedling traits at the three seedling densities. These results were consistent with the observation that all the traits displayed continuous variation in the RIL population, implying that seedling traits in rice were contributed by many loci with small effects. Furthermore, the existence of specific QTLs at different seeding densities also revealed the differences in seedling traits of indica rice and japonica rice in response to seeding densities.

Genome-wide additive and epistatic effects
Although the main-effect QTLs of young seedlings traits related to machine transplanting under different seeding densities had been fully explored. We did not have a complete understanding of the contributions of the '9311' and 'Nip' genomes to these traits at wholegenome level. Genome-wide additive and epistatic effects for each trait at three seeding densities were estimated based on the high-density bin map. The genome-wide additive and epistatic effects and the distribution of QTLs in Ex1 were detailed showed in Fig. 3. Trait. SH, seedling height (cm); FLSL, first leaf sheath length (cm); FLL, first leaf length (cm); SLL, second leaf length (cm); RDW, root dry weight (mg); SDW, shoot dry weight (mg). b Seeding density, LD, MD, and HD represented low density, medium density, and high density, respectively. c SD, standard deviation. d T-test between two parents, * and **, significant at P ≤ 0.05 and P ≤ 0.01. e Broad-sense heritability (%). f G × E, interaction of genotype and seeding density, * and **, significant at P ≤ 0.05 and P ≤ 0.01 Single-locus positive and negative additive effects were extensively distribution among all the seedling traits, however, the contributions of the allele from '9311' and 'NIP' differed across the traits and seeding densities (Additional file 1: Table S5). One of the most noticeable results was that the number of negative additive effect bins was more than that of positive additive effect bins for most traits, and the trend was consistent under different seeding densities (Table 3). However, there were a few exceptions, the number of positive and negative additive effect bins was approximately equal for SH and FLSL in Ex2 and Ex3. Besides, there were no bins with significant positive effect for some traits at both three seeding densities, such as FLL in all three experiments, SLL in Ex1, and RDW in Ex2 (Table 3). These results implied that the allele from '9311' contributing to increase seedling traits performance were widely genomic distribution. The comparisons of the distribution of additive effects under three seeding densities showed that the number of bins simultaneously detected under the three seeding densities took the advantage for all the other traits, except for RDW, the number of bins detected at the LD condition dominated. Another noteworthy result was that the bins that detected under HD were all contained under MD for FLL, while the bins detected under MD were contained under LD for RDW, and the bins detected under HD were included under both LD and MD for SLL (Fig. 4).
Digenetic interactions of each bin pairs across the entire genome at 0.0001 probability level were shown in Additional file 1: Table S6. Less numbers of digenic interactions contributing for seedling traits, explaining lower proportions of phenotypic variance than singlelocus analysis (Table 3). However, most of the interaction pairs individually accounted for more than 6% of the genotypic variation, which were not less than some main-effect QTLs. Significant interactions involve many marker loci, most of which were not detected in the single-locus analysis. However, there were still some QTLs found to interact with two or more other locus. The qFLSL 12.1 was identified simultaneously interacted with bin1368, bin1845 and bin2499 at HD. The qFLL 6.1 was detected to interacted with qFLL 2.2 and qFLL 9.1 at HD, and they were both showed positive epistasis effects. However, the interactions might change under different seeding densities. For example, the qFLL 6.1 was detected to interacted with bin1938 at LD in Ex1. Similarly, qSLL 2.1 at QTL hotspot bin411 interacted with bin2640 at LD, while it interacted with bin1290 at HD in Ex1. QTL for SDW at bin90 interacted with bin892 at LD and MD, while it interacted with bin431 and bin1302 at HD in Ex1 (Fig. 3). Another result worth noting was that the number of interactions with negative epistasis effects were predominant for SH, FLSL, SLL, RDW and SDW at three seeding densities, except for FLL (Table 3 and Additional file 1: Table S6). Taken together, genome-  wide additive and two-locus epistasis suggested a dynamic of the genetic control underlying different seeding densities.

Discussion
Our conclusions implicated the complex of genetic basis of rice seeding traits under different seeding densities which was only beginning to be acknowledged. In considering the results previous researches, we realized that the fundamental genetic mechanisms of seedling trait in rice were so complex, and partially caused by the diversity of rice germplasm [37,38,43]. Therefore, we compared identified QTL regions with previously reported QTL loci in other mapping populations and germplasm materials. qRDW 1.1 were located at bin72, which was previously reported affecting root number and root dry weight and length of mesocotyl [46][47][48]. The QTL hotspot at bin311 (qSH 1.1 , qFLSL 1.1 ) were overlapped with QTLs related to shoot length in previous reports and a candidate gene Os01g0904700 was identified responsible for increasing shoot length [35,37,49]. However, another main-effect QTL on chromosome 3 reported to increase seedling height and leaf sheath length were identified as OsGA20ox1, which was not detected in our research [41,49]. The QTL qSDW 5.1 were overlapping with the corresponding seedling vigor QTL mapped to the same locations [50,51]. Besides, another QTL hotspots on chromosome 6 (qFLSL 6.1 , qFLL 6.1 and qSLL 6.1 ) and qFLL 9.1 coincided with those previously reported QTLs related to shoot dry weight and shoot length in both directions of additive effects and the locations [37]. The co-location of QTLs for rice seedling traits indicated the reliability of these QTLs with common effects and could be found in different studies even in different populations or germplasm resources. One of the highlights of our research was conducting the seeding density experiments directly on the rice seedling tray, which was closely related to rice production. Although the quality of rice trays and substrates were strictly controlled in our research, the fluctuations of light, temperature and other climatic factors between years might still impact seedling growth. Nevertheless, some intriguing findings had been revealed. The growth of rice seedlings was significantly influenced by seeding density, as indicated in phenotypic data (Table 1). It showed that the seedling traits of the two parents exhibited similar trends as the seeding density increased, including SH, FLSL and SLL increases, SDW and RDW decreases. However, the differences of RDW gradually decreased, while SH, FLSL and SDW gradually increased as the seeding density increased, indicating the difference between aboveground and underground parts response to seeding density between '9311' and 'NIP'. The average traits performances of RIL population showed similar trends of the parents as the seeding density increased, and normally distributed on the seedling tray. Moreover, the interactions between seeding densities and genotypes had a greater influence on seedling traits indicating that the gene action related to seedling traits should be evaluated under certain densities. QTL analysis results indicated that '9311' contributed more QTL loci enhancing trait performances than 'NIP' for all six traits. The results of whole-genome additive effects under three seeding densities also supported that more  alleles responding to seeding density can be found in indica subspecies than in japonica subspecies. Although there were considerable overlaps in the genome-wide distribution of significant additive effects for the traits under the three seeding densities, such as SH, FLSL, FLL and SLL, special loci at certain seeding density still existed, which suggested that these loci were sensitive to seeding density and worthy further investigation (Fig. 4). The prevalence of negative epistasis effects indicated that the complementary two-locus homozygotes may not have marginal advantages over the means of the two parental genotypes. Moreover, some QTLs were found to interact with two or more loci, which suggested the main effects of QTLs were likely to be embedded in the interactions the rest of whole genome. Some interactions were specifically detected at certain seeding density, suggesting epistasis effect were also sensitive to seeding density. Thus, it may not be appropriate to interpret the single-locus marginal effects without specifying the genotypes and seeding density of the counterpart. Another highlight was that several QTLs were specifically detected at certain seeding density ( Table 2). Among 37 QTLs for six traits detected under the three seeding densities in our study, only 12 QTLs were detected in both three seeding densities. In other words, the rest QTLs were seedling density specific. This strongly suggested that QTL detection for seedling traits depended on the specific seeding density. For SH, all the three QTLs (qSH 1.1 , qSH 4.1 and qSDW 11.1 ) showed consistent effects on the traits across all three seeding densities, implying that the effect of seeding density on seedling height had universal effects. However, for FLSL, only  Fig. 3 The comparisons of genome-wide distributions of QTLs, additive and epistasis effects for seedling traits under three seeding densities in Ex1. a seedling height; b first leaf sheath length; c first leaf length; d second leaf length; e root dry weight; f shoot dry weight. The outermost circle represented the distribution of bin markers on 12 chromosomes. The histogram plots from the second to forth circle were the distributions of additive effects under LD, MD and HD, respectively. The blue, red and green bars represented the bins with significant additive effects under LD, MD and HD, respectively, and the black bars indicated non-significant additive effects. The blue squares, red circle and green triangle indicated the locations of QTLs detected under LD, MD and HD, respectively. The blue, red and green lines in the cycle indicated significant interactions under LD, MD and HD, respectively one QTL (qFLSL 1.1 ) was detected across the three seeding densities. qFLSL 5.1 was specifically detected at bin1281 region under MD and HD, while qFLSL 6.1 and qFLSL 12.1 specifically detected at LD and MD. Besides, one QTL hotspot for FLL and SLL at bin411 were detected at MD and HD, except for LD, and the allele of '9311' increased the trait values. Two QTLs for RDW, qRDW 1.1 and qRDW 3.1 were both detected under LD and MD, rather than HD. Among all the QTL identified for SDW (qSDW 2.1 , qSDW 4.1 , qSDW 5.1 , qSDW 5.2 and qSDW 11.1 ), qSDW 2.1 , qSDW 4.1 and qSDW 11.1 showed consistent effects on the traits across all three seeding densities, while qSDW 5.1 and qSDW 5.2 only detected at LD and MD, respectively. However, none of QTLs for RDW overlapped among the three seeding densities. Root dry weight were more easily affected by seeding density compared with shoot related traits, suggesting that there might be different mechanisms involved in root and shoot development under certain seeding density. However, the correspondences between those specifically detected QTLs and responses to the seeding density could not be clarified. Although the physiological and biochemical mechanisms of these QTLs have not been elucidated, specific QTLs might be useful for rice seedling improvement during mechanical transplanting.
Rice breeding strategies had been focused on yield, insects and diseases resistance, high nutrient efficiency and drought resistance in recent decades [52], less attention had been paid to improve the efficiency of mechanical transplanting through seedling traits improvements. The selection strategies of optimal seeding density for hybrid rice and conventional rice varieties were different from each other during mechanical transplanting [5]. Under high seeding density, indica rice grew faster than japonica rice seedlings, the decrease of seedling quality leads to increase of high-order tillering and yields decline [53], which supported our results that there were more QTLs and additive loci sensitive to seeding density in '9311' as compared with 'NIP'. Due to the higher seed cost of hybrid rice and seedling heterosis [54], the advantages of hybrids can only be exerted by low-density seeding, so the seedling quality should be improved as much as possible under the low-density seeding condition to ensure the completion of machine transplanting [17,26]. Conventional indica rice variety required highdensity seeding because of its lower seed cost and poor  SH  LD  36  488  12  42  80  99  0  33  166  363  4  50   MD  62  529  13  35  82  86  0  41  137  228  2  64   HD  82  465  6  42  73  109  1  53  166  159  3  79   FLSL  LD  104  225  5  12  168  294  1  79  173  183  1  88   MD  85  435  9  35  271  245  0  100  161  257  3  91   HD  74  382  12  19  351  243  0  118  138  267  2  120   FLL  LD  0  323  22  9  0  335  8  36  0  466  15  25   MD  0  497  16  9  0  278  25  17  0  553  18  31   HD  0  535  11  13  0  202  3  53  0  379  6  34   SLL  LD  0  437  4  29  39  321  1  40  21  215  4  62   MD  0  547  18  27  48  391  0  68  26  252  8  103   HD  0  491  10  31  19  461  0  68  23  381  9  99   RDW  LD  1  425  1  26  0  650  1  5  6  125  5  11   MD  0  345  2  4  0  184  0  3  0  99  2  9   HD  0  471  0  18  0  300  0  12  17  83  0  2   SDW  LD  24  484  8  41  63  334  6  15  74  282  8  29   MD  1  525  5  18  42  108  1  12  19  559  5  14   HD  5  416  2  55  21  97  4  14  35  440  7  18 a Trait, SH, seedling height; FLSL, first leaf sheath length; FLL, first leaf length; SLL, second leaf length; RDW, root dry weight (mg); SDW, shoot dry weight (mg). b Seeding density, LD, MD, and HD represented low density, medium density, and high density, respectively. c, d PAE and NAE, indicated the number of bins with positive additive effect and negative additive effect, respectively. Positive values indicate that alleles from Nip were in the direction of increasing the trait scores, and negative values indicate that alleles from 9311 are in the direction of increasing the score. e, f PEE, NEE indicated the number of two-locus interactions with significant positive and negative epistatic effects, respectively seedling quality than hybrid rice [5]. However, high density lead to further decrease of individual quality, the alleles and QTLs that exerted their functions under high-density seeding conditions were important for conventional indica rice seedling quality improvement. In the present study, QTLs specifically detected at low or high seeding densities may be important potential targets for improving seedling traits to adapt to machine transplanting. The interweaving force of rice seedlings was an important limiting factor for hybrid rice mechanical transplanting, which required enough root system to form firm structures on the premise of low seeding density [7]. Therefore, alleles and QTLs that functioned on root related traits under low seeding density were the best candidates to improve the seedling quality of hybrid rice seedlings. The qRDW 1.1 would be an ideal target for hybrid rice seedling improvement. The allele from '9311' strongly increased the root dry weight, as the seeding density decreased (Fig. 5a and b). Another potential target was qFLSL 5.1 , which was specifically detected at high seeding density, and the allele of '9311' increased the first leaf sheath length ( Fig. 5c and d). Although both the allele of 'NIP' and '9311' increased first leaf sheath length as seeding density increased, the '9311' allele was more sensitive to seeding density than the 'NIP' allele, resulting in a significant increase in FLSL of '9311' with the increase of seeding density. Thus, the allele of 'NIP' could improve first leaf sheath length for conventional indica rice varieties under high-density seeding conditions. Those QTLs detected under specific seeding density conditions may be involved in the shade avoidance response of rice seedlings. Detailed analysis of the regions of these QTLs revealed several candidate genes, involving in phytohormone signals, including abscisic acid (ABA), gibberellin, auxin, and ethylene and some abiotic-stresses related genes (Table 4). qRDW 1.1 were identified on chromosome 1. This region was found to contain five putative genes, two of them were involved in auxin (LOC_Os01g13520) [55,56] and ABA (LOC_ Os01g13530) [57] signals. NLP3 (LOC_Os01g13540) played a major role in nitrate uptake, translocation and QTLs specifically detected at certain seeding density. a lod curves and peak of qRDW 1.1 ; b Comparisons of root dry weight between 9311 and Nip genotypes under three seeding densities. c lod curves and peak at three seeding densities of qFLSL 5.1 ; d Comparisons of first leaf sheath length between 9311 and Nip genotypes at three densities. Ns, * and ** indicated no significant difference, significant at P ≤ 0.05 and P ≤ 0.01, respectively signaling in rice root [58]. Both LOC_Os01g13570 [59] and LOC_Os01g13740 [60] participated in abiotic stress response. Similarly, in the qFLL 2.1 region, three candidate genes including OsbZIP19 (LOC_Os02g14910), OsCNGC1 (LOC_Os02g15580) and ABA receptor PYL3 (LOC_Os02g15640) were identified. OsbZIP19 was upregulated in overexpression of ONAC022 plants and played a role against abiotic stresses in rice [61]. OsCNGC1 [62] and PYL3 [63] were both triggered by ABA signals. Four genes were found within the qFLSL 5.1 region. OsABI4 (LOC_Os05g28350) was the link between ABA and cytokinin signal transduction [64], likely to be the best candidate in this region. LOC_Os05g28500, a pentatricopeptide-repeat protein (PPR), specifically expressed in rice seedling and 20-day-old leaves, induced under the biotic and abiotic stresses [65]. LOC_ Os05g28730, encoding a C3HC4-type zinc finger domain containing protein [66] and LOC_Os05g28740 (universal stress protein) [67] were both reported to involved in abiotic stress response. Detailed functional analysis of these genomic regions and candidate genes through construction and evaluation of near-isogenic lines would further improve the understanding of genetic basis of seedling traits related to seeding densities as well as putting into practice according to variety types and cropping systems during mechanical transplanting.

Conclusions
In this study, we employed the technique of bowlshaped blanket-like seedling to assess seedling traits variations of a RIL population derived from '9311' and 'Nipponbare' under three seeding densities. Our results revealed the differences between indica rice and japonica rice seedling traits in response to seeding density. Several genomic regions containing QTL hotspots involved in different traits and specific QTLs (such as qRDW 1.1 and qFLSL 5.1 ) in diverse seeding densities had been detected. Genome-wide additive and two-locus epistasis suggested a dynamic of the genetic control underlying different seeding densities. Novel QTLs functioned under specific seeding density could be potential targets for marker-assisted selection (MAS) in rice breeding during mechanical transplanting.

Plant materials
The genetic population for QTL mapping in this study consisted of 213 recombinant inbred lines (RILs) derived by single-seed descent from a cross between Oryza sativa ssp. indica cv. 9311 and Oryza sativa ssp. japonica cv. Nipponbare (abbreviated to NIP) [68]. The original seeds come from the State Key Laboratory of Crop Genetic Improvement, National Plant Gene Research Center, Huazhong Agricultural University, Wuhan, China. All RILs and the parents were planted at Fuyang, Hangzhou in 2017. Mature seeds were harvested and stored at 4°C for use.

Experimental design and seedling raising
The experiment was carried out in Fuyang, Hangzhou, China Rice Research Institute. In this study, the seedling raising method was mechanical planting technique of bowl-shaped blanket-like seedling. The substrate used for rice seedlings raising is commercial substrate in the market without adding other ingredients. The seedlings were raised on plastic nursery tray with 58 cm long, 28 cm wide and 2.5 cm height, containing 420 holes. The seeding densities were set as three levels, 2, 4 and 6 grains per hole, and called low seeding density (LD), medium seeding density (MD) and high seeding density (HD), respectively. The experiments were repeated three times on May 2018, September 2018, and May 2019, referred to as Ex1, Ex2 and Ex3, respectively. The experiment was arranged as a split-plot design. The experiment was laid out in 3 blocks of 215 main plots, each split into 3 sub-plots. The RILs and two parents were applied to the main plots and the seeding density conditions to the sub-plots. According to the experiment design, plump seeds of each line were selected to detect germination rate, and then the qualified seeds were soaking with 25% prochloraz agent. After 2 days, the seeds were transferred to greenhouse incubator to accelerate germination for 12 h. The corresponding amounts of germination grains were measured for each seedling tray, and then evenly distributed on the surface of the nursery tray by the rice precision seeder, with equal amount of substrate on the back cover. Conventional field seedling management was adopted for seedling management.

Phenotyping of seedling traits
After 20 days of growth, the seedlings along with substrates for each genotype and seeding density in area of 10 cm × 10 cm in the middle of the nursery tray were selected and harvested. The sampling was repeated three times for each sub-plot. The phenotypes measured including seedling height (SH, cm), first leaf sheath length (FLSL, cm), first leaf length (FLL, cm), second leaf length (SLL, cm), root dry weight (RDW, mg) and shoot dry weight (SDW, mg). To determine RDW and SDW, the seedlings with substrate were dug out, placed in a sieve, and gently rinsed until no substrate remains. The average of the thirty seedlings with uniform size was treated as trait value in each sub-plot. The roots and shoots were then separated by cutting from the basal part of shoots. All the other traits were measured in the laboratory immediately following harvest, except shoot and root dry weight were measured after oven drying at 70°C for 3 days until constant weight.

Statistical analysis
A linear mixed model was used to estimate mean of each line, with lines as fixed effects, replications and blocks as random effects and in CropStat (v7.2.2007.3) (http://bbi.irri.org/products). The phenotypes and correlation coefficients among the traits were analyzed using R Statistics (R version 3.6.1) [69]. Analyses of variance (ANOVA) were used to test the interactions of genotype × environment (G × E) among three seedling densities. in each experiment. Broad-sense heritability (H) was calculated by an ANOVA using the formula: where σ 2 G , σ 2 GE and σ 2 ε were the genotypic variance, genotype × environment interaction variance and residual error variance, respectively, and r was the number of replicates and n was the number of experiments.
The single-locus additive effect of each bin across the whole genome for each experiment was calculated by the half of the difference between the means of the two homozygotes '9311' and 'NIP'. Any bin with heterozygous genotype was treated as a missing value. Analysis of variance with a threshold of F-value of 3.89 (P ≤ 0.05) were then used to test the significance of the additive effect at each bin. The bins with significant additive effects were then confirmed with 1000 permutation tests. Those bins were regarded as significant at P ≤ 0.05, unless no more than 5% of the random F-values were larger than the F-value from the original data [54]. Two-locus epistasis for each seeding density using all possible twolocus (bin) combinations in the RIL population were resolved using two-way ANOVA with a threshold of Fvalue of 11.13 (P ≤ 0.001). The calculation was based on unweighted cell means and the sums of squares were multiplied by the harmonic means of the cell sizes to form the test criteria [70]. Those significant digenic interactions were then confirmed by 10,000 random permutation tests, and the resulting 10,000 F-values were compared with the original Fvalue. If no more than one F-value from the random permutations was larger than the original F-value, the digenic interaction was regarded as significant (P ≤ 0.0001) [54,71,72]. Since the RIL populations were homozygous, the epistatic effects were characterized as additive by additive interaction.

QTL analysis
The complete linkage map of the RIL population constructed through high-density SNP marker analysis had been previously reported [68,73]. In our study, the RILs were genotyped based on SNPs generated from the whole-genome resequencing. The recombination maps of the RILs were aligned and compared for their genotypes for a 100-kb interval, resulting in a high-density  [73]. Each bin was then treated as a genetic marker for linkage map construction using MAPMAKER/EXP version 3.0b [74], resulting in a genetic linkage map of 1564.4 cM in length, covering all 12 chromosomes (Fig. 6 and Additional file 1: Table S1). QTL analysis was performed with the composite interval mapping (CIM) procedure of R/qtl function cim [75]. For the high-density bin map, because the bins were different from the traditional molecular markers, the scan window of R/qtl function cim was set to zero [76]. The likelihood ratio statistic was computed for each bin. Significant (P ≤ 0.05) LOD thresholds were determined by 1000 permutation tests. One QTL whose LOD value and the percentage of phenotypic variation explained exceeded 3.0 and 10% was declared as a main-effect QTL, and otherwise a minor QTL. The location of each QTL was determined based on the LOD peak location and a 1.5 LOD-drop support interval was calculated for each QTL to obtain a 95% confidence interval [76]. The additive effect and proportion of the phenotype variance explained by each QTL were determined by the linear model using the R [69].
Additional file 1: Table S1.Genotypes for all 2778 bins for the 213 RILs from the 9311/Nipponbare cross based on high quality SNPs. AA, Nipponbare genotype; BB, 9311 genotype; AB, heterozygous genotype. The position represents the physical position of each bin. Table S2. Phenotypic variation of the RIL population of 9311/Nipponbare cross among three seeding densities. SH, seedling height (cm); FLSL, first leaf sheath length (cm); FLL, first leaf length (cm); SLL, second leaf length (cm); SDW, shoot dry weight (mg); RDW, root dry weight (mg). Table S3. Correlation analyses among different seedling traits. LD, low density; MD, medium density; HD, high density; Ex1, Ex2 and Ex3 indicated the three experiments, respectively. SH, seedling height; FLSL, first leaf sheath length; FLL, first leaf length; SLL, second leaf length; SDW, shoot dry weight; RDW, root dry weight. * and ** indicated significant at P ≤ 0.05 and P ≤ 0.01, respectively. Table S4. The QTLs identified for seedling traits from the RIL population of 9311/Nipponbare cross under three seeding densities. LD, low density; MD, medium density; HD, high density; SH, seedling height; FLSL, first leaf sheath length; FLL, first leaf length; SLL, second leaf length; SDW, shoot dry weight; RDW, root dry weight. The position is the lod peak of each QTL. The interval is 1.5-LOD support interval of the QTL. Positive and negative additive effect values indicate that the allele from NIP and 9311 increase trait values, respectively. The Var is variation (%) explained by each QTL. Table S5. Additive effects of each bin under three seeding densities. A_score, additive effect, in which positive values indicate that alleles from NIP are in the direction of increasing the trait scores, and negative values indicate that alleles from 9311 are in the direction of increasing the score. NA, NIP and 9311 indicated the additive effect non-significant, positive significant and negative significant, respectively. LD, low density; MD, medium density; HD, high density; SH, seedling height; FLSL, first leaf sheath length; FLL, first leaf length; SLL, second leaf length; SDW, shoot dry weight; RDW, root dry weight. Table S6. Two-locus combinations showing significant epistatic effects under three seeding densities. LD, low density; MD, medium density; HD, high density; SH, seedling height; FLSL, first leaf sheath length; FLL, first leaf length; SLL, second leaf length; SDW, shoot dry weight; RDW, root dry weight. Var, the percentage (%) of variation explained by the interaction.   Availability of data and materials All data generated or analysed during this study are included in this published article and its supplementary information files.
Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.