{'fraction_non_alphanumeric': 0.04417928222980556, 'fraction_numerical': 0.029384146760414998, 'mean_word_length': 4.315802288775262, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 9, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A combinatorial strategy to identify various types of QTLs for quantitative traits using extreme phenotype individuals in an F2 population

Theoretical and applied studies demonstrate the difficulty of detecting extremely over-dominant and small-effect genes for quantitative traits via bulked segregant analysis (BSA) in an F2 population. To address this issue, we proposed an integrated strategy for mapping various types of quantitative trait loci (QTLs) for quantitative traits via a combination of BSA and whole-genome sequencing. In this strategy, the numbers of read counts of marker alleles in two extreme pools were used to predict the numbers of read counts of marker genotypes. These observed and predicted numbers were used to construct a new statistic, Gw, for detecting quantitative trait genes (QTGs), and the method was named dQTG-seq1. This method was significantly better than existing BSA methods. If the goal was to identify extremely over-dominant and small-effect genes, another reserved DNA/RNA sample from each extreme phenotype F2 plant was sequenced, and the observed numbers of marker alleles and genotypes were used to calculate Gw to detect QTGs; this method was named dQTG-seq2. In simulated and real rice dataset analyses, dQTG-seq2 could identify many more extremely over-dominant and small-effect genes than BSA and QTL mapping methods. dQTG-seq2 may be extended to other heterogeneous mapping populations. The significance threshold of Gw in this study was determined by permutation experiments. In addition, a handbook for the R software dQTG.seq, which is available at https://cran.r-project.org/web/packages/dQTG.seq/index.html, has been provided in the supplemental materials for the users’ convenience. This study provides a new strategy for identifying all types of QTLs for quantitative traits in an F2 population.


INTRODUCTION
In the past several decades, many methods have been established to mine elite genes for quantitative traits in animals, plants, and human beings; these include quantitative trait locus (QTL) mapping (Lander and Botstein, 1989), genome-wide association studies (Zhang et al., 2005;Yu et al., 2006), and bulked segregant analysis (BSA) (Giovannoni et al., 1991;Michelmore et al., 1991;Zhang et al., 1994). In real data analysis, it is difficult to identify small-effect and linked genes (Kroymann and Mitchell-Olds, 2005;Mackay et al., 2009;Wang et al., 2016;Wen et al., 2019). In theory, it is also difficult to detect extremely over-dominant genes in BSA (Takagi et al., 2013a;Schneeberger, 2014;Haase et al., 2015;Wang et al., 2019). As we know, quantitative traits are controlled by both a few major genes and a series of polygenes, whereas over-dominant genes can be used to explain heterosis, which is a natural phenomenon. Thus, it is necessary to investigate methods for detecting extremely overdominant and small-effect genes in BSA (Cockram and Mackay, 2018). developed, and they have been summarized in Li and Xu, 2021. These methods are mainly used to detect genes for qualitative traits, mutants, and quantitative traits. In the early stage, BSA methods were widely used to identify dominant genes for qualitative traits and then increase molecular marker density in each F 2 individual to fine-map target genes such as the dominant white locus in chicken (Ruyter-Spira et al., 1997). Recently, microsatellite markers were replaced by single nucleotide polymorphism (SNP) markers in the references of Li et al. (2018), Feng et al. (2019), and Liu et al. (2019). BSA methods were then widely used to identify mutant genes for qualitative traits. Schneeberger et al. (2009) first integrated BSA with whole-genome resequencing in a large pool of mutant F 2 plants to rapidly identify recessive/dominant mutant genes based on the frequency of mutant alleles using SHOREmap software. When a smaller number of mutant plants was available, Cuperus et al. (2010) used wild-type SNP enrichment in a pool of mutant F 2 plants to identify mutant genes. Austin et al. (2011) used the frequency of wild-type SNP alleles in a pool of mutant F 2 plants to locate recessive mutant genes by the next generation mapping method, and Lindner et al. (2012) used the causative SNP/non-causative SNP ratio in a pool of mutant BC 2 individuals to clone lethal mutations. For species with large and incomplete reference genomes, CloudMap (Minevich et al., 2012), SNPtrack (Leshchiner et al., 2012), Fishyskeleton (Bowen et al., 2012), MegaMapper (Obholzer et al., 2012), and MMAPPR (Hill et al., 2013) are available. When a mutant is dominant, the reverse mapping and F 3 screen strategies of Smith et al. (2016) are available. Abe et al. (2012), Fekih et al. (2013), and Takagi et al. (2013b) used the SNP index in a pool of mutant F 2 (M 3 ) plants from a cross between a mutant and its wild type to identify mutant genes with MutMap, MutMap+, and MutMap-Gap, respectively. To simultaneously identify multiple causal mutations, SIMM (Yan et al., 2016) and BSA-seq (Wang et al., 2021) are available. However, these recessive/dominant mutant traits are different from the dominant effects of quantitative traits.
Finally, BSA methods are widely used to identify quantitative trait genes (QTGs). As we know, quantitative traits are controlled by a few major genes and a series of polygenes with small genetic effects, making polygene identification difficult (Kroymann and Mitchell-Olds, 2005;Mackay et al., 2009). To identify these polygenes, SHOREmap (Schneeberger et al., 2009), MMAPPR (Hill et al., 2013), MutMap (Abe et al., 2012), MutMap+ (Fekih et al., 2013), and MutMap-Gap (Takagi et al., 2013b) are available. Ehrenreich et al. (2010) used allelic frequency to identify yeast polygenes via BSA-seq. The Z-statistic is also used to detect allelic frequency differences at the target locus in high and low pools (Huang et al., 2012;Benowitz et al., 2019). If several markers in one window (or region) are used to calculate the average Z-statistic, this is the Z 0 -statistic (Haase et al., 2015). Similarly, a smoothed G 0 of the standard G statistic was proposed in Magwene et al. (2011). In addition, a dynamic Bayesian network via MULTI POOL (Edwards and Gifford, 2012) is used to capture the association among nearby loci in a pooled sequencing experiment; a hidden Markov model (HMM) is used to calculate the most probable states via EXPLoRA, indicating genomic regions that may be likely to contain trait-related genes (Duitama et al., 2014;Claesen and Burzykowski, 2015); a nonhomogeneous HMM with a transition matrix is used to detect the association of markers with genes (Ghavidel et al., 2015). All the approaches above make use of high and low pools. To increase the power of QTG detection, Wang et al. (2019) added an additional, middle pool in GradedPool-Seq to identify the rice heterotic QTL GW3p6. If no secondary F 2 populations are available, a primary F 2 population from a cross between any two homozygous parents may be used. In this case, the DSNP index is also used to detect QTLs, and the procedure is termed QTL-seq (Takagi et al., 2013a). Clearly, the results may be affected by genetic background, although the method may be effective in some instances. To address this issue, Zhang et al. (2019) proposed a new design and a smooth-LOD (logarithm of odds) statistic. In analyses of real data, however, it is difficult to detect extremely overdominant and small-effect genes for quantitative traits (Ehrenreich et al., 2010;Takagi et al., 2013a;Schneeberger, 2014;Haase et al., 2015;Wang et al., 2019).
A dominant mutant of a single-gene qualitative trait is different from the dominant effect of a QTL. In theory, the former is a qualitative trait, and the latter is a quantitative trait; they have different genetic structures in the two extreme phenotype pools. The structure of high and low pools is QQ + Qq via qq for a dominant mutant and QQ + qq via Qq for a fully over-dominant (zero additive and non-zero dominance) QTL for a quantitative trait (Haase et al., 2015). Thus, the allelic frequency of a homozygous genotype pool is indistinguishable from the allelic frequency of a heterozygous genotype pool in the fully over-dominant QTL model, and the fully over-dominant QTL becomes undetectable (Supplemental File 1). In the same way, it is difficult to detect extremely over-dominant QTLs (Schneeberger, 2014;Wang et al., 2019).
To address the issues above, we propose a combinatorial strategy for identifying various types of QTLs for quantitative traits in the F 2 population by a combination of BSA and whole-genome sequencing. Specifically, equally mixed DNA sample of all plants in each extreme pool was sequenced. On the basis of the observed numbers of read counts of marker alleles in the two extreme pools, we propose a method for predicting the numbers of read counts of marker genotypes. The observed allele and predicted genotype numbers were used to construct a new statistic, G w, to identify marker-trait associations. If the goal was to detect extremely over-dominant and small-effect genes, another reserved DNA/RNA sample from each extreme phenotype individual was sequenced one plant at a time, and the observed numbers of marker alleles and genotypes were used to calculate G w to identify marker-trait associations. The former strategy is a new BSA method named dQTG-seq1, and the latter is a new individual segregant analysis named dQTG-seq2, where ''d'' indicates the dominant effect of the QTG (see materials and methods). Real and simulated datasets were used to analyze the sensitivity and specificity of the new strategy, and we also investigated the best sampling plan for the new strategy.

RESULTS
A combinatorial strategy for mapping all types of QTLs for quantitative traits in an F 2 population by a combination of BSA and whole-genome sequencing The new strategy in this study was proposed to quickly identify QTGs for quantitative traits in an F 2 population, especially for extremely over-dominant and small-effect genes ( Figure 1). There were six types of F 2 designs available. They included an F 2 from a cross between any two homozygous parents (F 2 -I), an immortalized F 2 (IMF 2 ) derived from a DH (doubled haploid) or RIL (recombinant inbred line) population (F 2 -II) (Gardiner et al., 1993;Hua et al., 2003), a secondary F 2 from a cross between a mutant and its wild type (F 2 -III) (Abe et al., 2012), a secondary F 2 from a cross between two near-isogenic lines (F 2 -IV), a secondary F 2 from the selfing of residual heterozygous lines of a targeted QTL (F 2 -V) (Yamanaka et al., 2005), and the improved secondary F 2 of Zhang et al. (2019) from the selfing of residual heterozygous lines of a targeted QTL (F 2 -VI). When one of these F 2 populations was constructed, the new strategy was available.
F 2 -I to F 2 -V have been described in previous studies; F 2 -VI here is slightly different from that in Zhang et al. (2019). The improved F 2 design involves the following steps: (1) Develop F 1 and F 2 populations from the cross between any two inbred lines. (2) Identify QTL for the quantitative trait of interest using traitphenotype and marker-genotype datasets from the F 2 population. (3) If the additive (a) and dominant (d) effects of the target QTL have the same sign, P 2 should be crossed with F 1 to construct BC 1(P2) F 1 . If they have opposite signs, P 1 should be crossed with F 1 to construct BC 1(P1) F 1 . This is because the large difference in the two genotypic values in BC 1(Ps) F 1 (s = 1, 2) can increase the power of QTL detection in BSA. To select BC 1(Ps) F 1 individuals with a heterozygous genotype at the target QTL and homozygous genotypes at non-target QTLs, molecular markers tightly linked to each detected QTL in the F 2 population are available. The selected BC 1(Ps) F 1 individuals are self-pollinated to produce a secondary F 2 population . To ensure a large population, multiple target BC s F 1 individuals should be selected. (4) The above-mentioned secondary F 2 individuals are measured for the trait of interest, and the plants with minimum and maximum phenotypes are selected separately within each family. This is the F 2 -VI design used in the present study.
In each aforementioned F 2 design, an equally mixed DNA/RNA sample of all individuals in each extreme pool is separately subjected to whole-genome sequencing with deep genome coverage c. The sequenced datasets for each mixed sample are mapped to the reference genome, and genomic variants across the genome are identified. On the basis of the observed numbers of read counts of marker alleles Q (j = 1) and q (j = 2) in low (i = 1) and high (i = 2) pools, we propose a method to predict the numbers of read counts of marker genotypes QQ (j = 1), Qq (j = 2), and qq (j = 3) in the two extreme pools (see materials and methods). Thus, the observed allele and predicted genotype numbers are used to calculate the new statistic G w where G 1 and G 2 are the standard G statistics of Magwene read numbers n Qij and their expectations E(n Qij ), and G 2 = P 2 i = 1 P 3 j = 1 n Gij ln½n Gij =Eðn Gij Þ using genotypic read numbers n Gij and their expectations E(n Gij ). The G w statistic was used to identify QTLs for quantitative traits, and the strategy was named dQTG-seq1.
Although dQTG-seq1 has higher power for QTL detection than existing BSA methods, it may fail to detect extremely and fully over-dominant and small-effect genes, for which all existing BSA methods including dQTG-seq1 do not work well. To address this issue, additional reserved DNA/RNA sample from each extreme phenotype individual is separately subjected to wholegenome sequencing one plant at a time. In this case, the observed numbers of marker alleles and genotypes in the two pools can be used to calculate G w to identify QTGs in a strategy we have named dQTG-seq2.
Monte Carlo simulation studies dQTG-seq1 has higher power than other BSA methods To demonstrate that dQTG-seq1 has higher power than other BSA methods, one small-effect (5%) QTL was simulated under four kinds of genetic models: additive (a s 0 and d = 0), completely dominant (a = d), and over-dominant (d = 2a and d = 3a) (Supplemental Table 1). All samples in simulated dataset I were analyzed by the dQTG-seq1, ED (Euclidean distance), G 0 , and DSNP methods. dQTG-seq1 had significantly higher power than the other BSA methods; its power was 10.0%, 8.5%, 4.7%, and 4.1% higher than those of the other BSA methods in situations with dominance ratios of 0.0, 1.0, 2.0, and 3.0, respectively ( Figure 2A). dQTG-seq2 can identify more loci with large jd =aj and small jDSNP indexj To investigate whether the dQTG-seq1 and dQTG-seq2 methods can identify extremely and fully over-dominant and small-effect genes in the F 2 population, one small-effect (5%) QTL was simulated under extremely (d/a = 4.0) and fully (a = 0 and d s 0) over-dominant models in a secondary F 2 (F 2 -III to F 2 -VI) population. Each simulation sample in simulated dataset II was analyzed by the dQTG-seq1, dQTG-seq2, ED, G 0 , DSNP, and ICIM methods, and the power, false positive rate (FPR), false negative rate (FNR), and accuracy of QTL position estimates in QTL detection were calculated. Under the extremely over-dominant model, the powers of QTL detection using dQTG-seq2 (71.3%) and ICIM (70.0%) were significantly higher than those using the dQTG-seq1 (14.3%), ED (12.7%), G 0 (13.7%), and DSNP (13.3%) methods. Similar trends were observed in the FPR, FNR, and accuracy of QTL position estimates ( Figure 2B; Supplemental Table 2). Under the fully over-dominant model, the powers of QTL detection using dQTG-seq2 (74.8%) and ICIM (73.0%) increased slightly, whereas the powers of the other methods decreased significantly (<5%). The trends in FPR, FNR, and accuracy of QTL position estimates were similar to those in the extremely overdominant model ( Figure 2B; Supplemental Table 2).
To assess the applicability of dQTG-seq1 and dQTG-seq2 to a primary F 2 population from a cross between any two homozygous parents (F 2 -I and F 2 -II), five QTLs were simulated under extremely (d/a = 4) and fully over-dominant (a = 0 and d s 0) models in simulated dataset III, and their sizes were set to 2.5%, 5.0%, 5.0%, 10.0%, and 10.0%. Under the extremely over-dominant model, the average powers of dQTG-seq2 (49.0%) and ICIM (30.1%) were significantly higher than those of the dQTG-seq1 (15.8%), ED (9.3%), G 0 (13.9%), and DSNP (13.3%) methods, and the FPR, FNR, and accuracy of QTL position estimates from dQTG-seq2 and ICIM were better than those from the other methods ( Figure 2C; Supplemental Table 3). Under the fully over-dominant model, the average powers of dQTG-seq2 (50.2%) and ICIM (28.4%) showed slight changes, whereas the average powers of the dQTG-seq1 (7.7%), ED (4.0%), G 0 (6.8%), and DSNP (6.6%) methods decreased. The trends in FPR, FNR, and accuracy of QTL position estimates were similar to those in the extremely over-dominant model ( Figure 2D; Supplemental Table 3). Interestingly, we found that QTL 4 and QTL 5 with a size of 10% had lower powers than QTL 2 and QTL 3 with a size of 5%. This may be because the effects of QTL 4 and QTL 5 were in opposite directions.
In conclusion, from the second and third simulation experiments, dQTG-seq2 can detect extremely and fully over-dominant and small-effect genes, but the other BSA methods, such as dQTG-seq1, ED, G 0 , and DSNP, do not work well.
To determine whether dQTG-seq1 and dQTGseq2 can be used to detect multiple QTLs on a chromosome, two 10% QTLs on one chromosome were simulated in an F 2 population; the genetic model was d = 0.5a, and all the samples in simulated dataset IV were analyzed by dQTG-seq1 and dQTG-seq2. The two QTLs could be detected (Supplemental Table 4), and the new methods can therefore be used to detect multiple QTLs on a chromosome.
Analyses of real data from maize and rice dQTG-seq1, as well as the G 0 and DSNP methods, can identify a QTG for maize plant height in a secondary F 2 population To validate dQTG-seq1 in a secondary F 2 population (F 2 -III to F 2 -VI), the maize dataset from Zhang et al. (2019) was reanalyzed. In Zhang et al. (2019), the numbers of read counts of marker alleles obtained from whole-genome sequencing were used to detect a QTN for plant height. The target QTL on chromosome 7 in the F 2 population was also detected in a secondary F 2 population, and its candidate gene around this QTN was confirmed by a molecular biology experiment . In the present study, we predicted the numbers of read counts of marker genotypes, and the observed allele and predicted genotype numbers were used to detect the QTN. In Figure 3, the QTN in Zhang et al. (2019) was also found in the secondary F 2 population by the dQTG-seq1, G 0 , and DSNP methods, indicating the feasibility of dQTG-seq1 for detection of target QTN in a secondary F 2 population. dQTG-seq1 identifies more known genes for rice yield in a primary F 2 population than the previous BSA (DSNP and G 0 ) methods To validate dQTG-seq1 in a primary IMF 2 population (F 2 -I and F 2 -II) from a cross between any two homozygous parents, the In Zhou et al. (2012), 11 QTLs were identified as significantly associated with rice yield. Around these QTLs, there were eight previously reported genes (Supplemental Table 5). In this study, 20% extremely low-yield plants and 20% extremely high-yield plants were grouped into low and high pools. In each pool, we calculated the numbers of marker alleles and predicted the numbers of marker genotypes. These numbers were used to calculate G w . As a result, 18 QTLs were found to be significantly associated with rice yield (Supplemental Table 5). Around these QTLs, 11 previously reported genes were found to be truly associated with rice yield (Supplemental Table 6; Figure 4A). When the previous (G 0 and DSNP) BSA methods were used, only eight previously reported genes were found by each method (Figure 4B and 4C). Clearly, dQTG-seq1 identified more known genes than the G 0 and DSNP methods, demonstrating the advantage of dQTG-seq1 over the previous BSA (G 0 and DSNP) methods. dQTG-seq1 and existing BSA methods identify fewer large jd =aj and small-effect known genes for rice yield in a primary F 2 population than CIM and ICIM To compare dQTG-seq1 with previous QTL mapping methods, all individuals in Zhou et al. (2012) were re-analyzed by the CIM and ICIM methods. As a result, eight (CIM) and six (ICIM) known genes were found to be truly associated with rice yield ( Figure 4D and 4E). Among these known genes identified by QTL mapping methods, five were not identified by the BSA methods. Here, we calculated the jd =aj and jDSNP indexj of the five known genes and found that four had a large jd =aj (3.04, 5.93, 4.22, and -6.46), and one had a small jDSNP indexj. This result highlights the difficulty of using the dQTG-seq1, G 0 , and DSNP methods to detect genes with jd =aj > 2 or small jDSNP indexj (%0.15). To confirm this viewpoint, we calculated the jd =aj and jDSNP indexj of the 11 known genes identified by dQTG-seq1. The value of jd =aj ranged from 0.00 to 1.73, and its average was 0.72 ± 0.55; jDSNP indexj ranged from 0.13 to 0.22, and its average was 0.19 ± 0.026 (Supplemental Table 6). The low jd =aj and large jDSNP indexj phenomena were also observed in the G 0 and DSNP methods. Thus, dQTG-seq1 and existing BSA methods have shortcomings in the detection of large jd =aj and small jDSNP indexj genes, although dQTG-seq1 has higher power than existing BSA methods.
Here, we provide some simulation and theoretical analyses to support the thresholds of jd =aj > 2 and jDSNP indexj % 0.15. First, in Figure 2A, jd =aj = 2.0 was viewed as a critical value based on the stability of power changes. That is, the power of small-effect (5%) QTL detection in 500 F 2 individuals is approximately 30%. Second, we carried out one theoretical analysis to explain why the critical value of jDSNP indexj is 0.15. In situations where QTL size ranged from 1% to 5% and sampling fraction was 5%, 10%, 15%, and 20% in the F 2 population, we calculated allelic frequency differences between the two extreme pools. We found that most allelic frequency differences were less than 0.15 (Supplemental Figure 1), in which it was difficult to detect these loci. The above two critical values were consistent with those in analysis of real data (Supplemental Table 6). dQTG-seq2 identifies more large jd =aj and small jDSNP indexj loci/genes for rice yield in an F 2 population To address the aforementioned issue with dQTG-seq1 and existing BSA methods, each extreme plant from the low and high pools of the primary or secondary F 2 population was sequenced, and the observed numbers of marker alleles and genotypes were used to calculate G w (Figure 1). To validate dQTG-seq2, the above-mentioned rice dataset from Zhou et al. (2012) was reanalyzed. All the IMF 2 individuals were sorted according to their phenotypic values, and 20% of the extreme high and extreme low individuals were selected to form the high and low pools. From the 20% extreme low and high pools, we could obtain 5%, 10%, and 15% extreme low and high pools. In these four datasets, dQTG-seq2 was used to detect QTLs for rice yield. As a result, 33 QTLs were found to be significantly associated with rice yield ( Figure 4F; Supplemental Table 5). Among these QTLs, 26 previously reported genes were found to be truly associated with rice yield, jd =aj ranged from 0.04 to 36.61 with an average of 4.37 ± 7.69, and there were eight known genes with jd =aj > 2.0, four and three of which were identified by CIM and ICIM, respectively. The jDSNP indexj ranged from 0.00 to 0.22 with an average of 0.13 ± 0.064, and there were 14 known genes with jDSNP indexj < 0.15, one of which was identified by CIM (Supplemental Table 6). These results confirmed the advantages of dQTG-seq2 over dQTG-seq1 and existing BSA methods.
Among the above-mentioned 27 known genes, 25 were validated by one-way ANOVA by calculating the average values of individuals with each genotype (Supplemental Figure 2).

Suitable sampling plans via Monte Carlo simulation studies
To obtain suitable sampling plans for the dQTG-seq1 and dQTG-seq2 methods, we conducted a series of Monte Carlo simulation studies. In these studies, sample sizes were set to 250, 500, 1 000, and 2 000; QTL sizes (r 2 , %) were set to 2%, 5%, and 10%; 30, 50, and 80 extreme individuals were sampled in each extreme pool. There were 1 000 replicates. All samples in simulated dataset V were analyzed by dQTG-seq1 and dQTG-seq2, and the results are shown in Figure 5 and Supplemental Table 7. If sampling plans in which the power of QTL detection is >40% were suggested to be suitable to users, some suitable sampling plans were obtained. That is, 80 extreme individuals in 2 000 F 2 plants were selected for each pool to identify 2% Horizontal dotted lines indicate thresholds for significant QTLs, where the critical values were determined by permutation experiments for the new methods and by the R package QTLseqr for the G 0 and DSNP methods and set to LOD = 2.5 for the CIM and ICIM methods. Various statistics of genomewide scanning using new and existing methods are indicated by blue line. The genes with jd =aj < 2.0, small effects, and jd =aj R 2.0 are indicated by black, red, and pink, respectively. If the dominance ratio of a small-effect gene is larger than 2.0, the gene name is in pink, and its corresponding solid line is in red. The dataset was derived from Zhou et al. (2012).
QTLs via dQTG-seq2, whereas more than 2 000 F 2 individuals were needed for dQTG-seq1. Thirty extreme individuals in at least 250 F 2 individuals were selected for each pool to identify 5%-10% QTLs via dQTG-seq1 and dQTG-seq2. In summary, selection of 80 extreme individuals for each pool from a large F 2 population was suitable for identifying small-effect (2%) genes, whereas selection of 30 extreme individuals for each pool from a general F 2 population was suitable for identifying large-effect genes.

DISCUSSION
Significant progress in BSA has been made in this study. First, a new BSA strategy was proposed to quickly detect QTLs for

Plant Communications
quantitative traits in primary or secondary F 2 population. Specifically, equally mixed DNA/RNA samples of all the plants in each extreme pool are sequenced one pool at a time. The sequencing results are used to identify QTLs via dQTG-seq1. If the goal is to identify large jd =aj and/or small jDSNP indexj genes, reserved DNA/RNA sample from each extreme plant is sequenced one plant at a time. The observed numbers of marker alleles and genotypes are used to detect QTLs via dQTG-seq2. In analyses of simulated and real data, dQTG-seq2 can identify extremely and fully over-dominant and small-effect QTLs, whereas dQTG-seq1 is feasible for detecting low jd =aj and large-effect loci. The two methods had significantly higher powers of QTL detection than existing BSA methods. This is because the new methods utilize both allelic and genotypic frequency differences between high and low pools ( Table 1). As shown in Figure 6A-6E, jDSNP indexj decreases as jd =aj increases, especially in the fully over-dominant model in which jDSNP indexj = 0, but differences in marker genotypic frequencies exist, demonstrating the effectiveness of detecting a gene or locus using dQTG-seq2. Then, in BSA, we predicted the numbers of read counts of marker genotypes from the numbers of read counts of marker alleles, and the observed allele and predicted genotype numbers were used to construct the new statistic G w . We found that the predicted values were acceptable when jd =aj < 2.0 and jDSNP indexj > 0.15, even in a primary F 2 population (Supplemental Tables 1 and 5; Supplemental Figure 1). If jd =aj R 2.0, jDSNP indexj < 0.15, a = 0, and d s 0, the predicted values were inaccurate, the power of QTL detection was low, and the FPR was high (>5W) (Supplemental Tables 1 and 5; Supplemental Figure 1), indicating the difficulty of detecting extremely and fully over-dominant and small-effect genes or loci via dQTG-seq1. We therefore proposed sequencing each extreme individual in the two pools ( Figure 1) so that the numbers of marker alleles and genotypes could be obtained, and extremely and fully over-dominant and small-effect genes could be detected by dQTG-seq2 with high power and accuracy (Figures 2 and 4; Supplemental Table 6) and a low false-positive rate (Figure 2; Supplemental Table 2). In the analysis of real data for rice yield, dQTG-seq2 identified eight over-dominant known genes, four of which were detected by CIM and ICIM; the other BSA methods found no known genes. More importantly, dQTG-seq2 identified eight small-effect known genes, whereas the other methods detected no known small-effect genes. Finally, dQTG-seq2 had some slight advantages over CIM and ICIM in simulation and real data analysis, such as the low cost of marker genotyping and slightly higher power and accuracy. As shown in Figure 4, dQTG-seq2 identified more known genes than CIM and ICIM, suggesting that CIM and ICIM can be replaced by dQTG-seq2 in QTL mapping.
In quantitative genetics, over-dominance and complete dominance are defined as jd =aj > 1.2 and 0.8 < jd =aj < 1.2, respectively (Stuber et al., 1987). In this study, extreme overdominance was defined as jd =aj R 4.0, and full overdominance was defined as d s 0 and a = 0. Small-and large-effect QTLs (genes) were distinguished between r 2 < 5% and r 2 R 5% based on the results in Supplemental Figure 1 and the analysis of real rice data. The definitions in this study are helpful for understanding the new methods.
The new methods can be used to identify QTLs in the aforementioned six types of F 2 population. As we know, a secondary F 2 population is better than a primary F 2 population for fine-mapping of QTLs. This is also true for BSA. If no QTLs in a secondary F 2 (F 2 -III to F 2 -VI) or no ideal QTLs in a primary F 2 (F 2 -I and F 2 -II) are detected by dQTG-seq1, we recommend sequencing another reserved DNA/RNA sample from each extreme plant, and the observed numbers of marker alleles and genotypes can be used to identify QTLs via dQTG-seq2. We should point out that F 2 -VI is an improved F 2 design based on Zhang et al. (2019), and the designs differ in whether two parents may be used to backcross the F 1 . In addition, dQTG-seq2 is theoretically applicable to any heterogeneous mapping population, such as BC 1 F s , F s (s R 3), genetic mating design, or association mapping population, because the numbers of marker alleles and genotypes are observed, and their theoretical numbers are calculated from a contingency table. Thus, dQTG-seq2 is feasible. Detailed results will be reported soon.
Sometimes, F 2 individuals have a large experimental error. First, the F 2:3 design (Zhang and Xu, 2004) can be used to overcome this problem. That is, DNA/RNA is extracted from each F 2 plant, Table 1. Mathematical differences among dQTG-seq1, dQTG-seq2, and other BSA methods.

Statistic Observed alleles Predicted number of genotypes Observed number of genotypes
, n Qij are the number of read counts of marker alleles, and E(n Qij ) is the expectation of n Qij ; G 2 = ; n Gij are the number of read counts of marker genotypes, and E(n Gij ) is the expectation of n Gij . P AL and P AH (P aL and P aH ) are the frequencies of marker allele A (a) in extreme low and high pools, respectively. k j is the weight of the jth marker in the window. and extreme individuals are determined based on the average phenotypes of F 2:3 families. This design is frequently used in maize and cotton genetic analysis. IMF 2 derived from DH lines or RILs is also used to overcome this problem (Hua et al., 2003). This is because replicated experiments can be implemented in IMF 2 , especially from RILs with more recombinants.
The new methods are different from previous BSA and QTL mapping methods. To reduce the time required to construct the mapping population and the cost of marker genotyping (Darvasi and Soller, 1994), BSA is proposed at an early stage to detect the association of a marker with a trait (Giovannoni et al., 1991;Michelmore et al., 1991). This method is also frequently used to locate a single-gene trait in Mendelian genetics (Ruyter-Spira et al., 1997). Specifically, some closely linked markers around one marker with significant differences in allele frequency between the two extreme pools are used to fine-map the gene via linkage analysis. With the development of sequencing technology and reductions in sequencing costs, the integration of BSA with sequencing technology provides new opportunities for quickly detecting QTGs. To date, a series of methods have been proposed (Magwene et al., 2011;Abe et al., 2012;Takagi et al., 2013aTakagi et al., , 2013bYan et al., 2016). In these methods, the numbers of read counts of marker alleles in low and high pools in primary or secondary F 2 populations are used to construct their statistics. In this study, the numbers of (read counts of) marker alleles and genotypes in the two extreme pools are used to construct the new statistic G w . If the trait is controlled by a fully over-dominant gene, in theory there are no differences in allelic frequencies of markers around the trait gene between the two extreme pools (Supplemental File 1), highlighting the inability of existing BSA methods to locate fully over-dominant genes (Schneeberger, 2014;Wang et al., 2019). Clearly, the present study solves this problem, which has been pending for some time.
As described in model 1, G w consists of G 1 and G 2 , which are based on the numbers of read counts of marker alleles and genotypes, respectively. If G 1 = 0, G 2 is available; if G 2 = 0, G w is G 0 of Magwene et al. (2011). Thus, G w is an extension of G 0 . Although the distributions of G 1 and G 2 are c 2 (Magwene et al., 2011), G 1 and G 2 are correlated; thus, the distribution of G w is unknown. This makes it difficult to determine the critical value of G w for significant QTGs. In the present study, 1 000 permutation experiments were used to determine the critical value. The results are presented in Supplemental Figure 3. They show that the critical value of G w is 6.52-7.24 at the 0.10 probability level and 7.60-8.52 at the 0.05 probability level using the dQTG-seq2 method. We suggest that the 0.10 threshold be used with dQTG-seq1 and the 0.05 threshold be used with dQTG-seq2. This is because the predicted numbers of read counts of marker genotypes have residual error. This may explain why different critical values in Figures 4A (6.52) and 4F (7.60) are used. In addition, we have listed additional G w thresholds under various probability levels in Supplemental Table 8 and Supplemental Figure 4.
The new method is also applicable in other situations. If structural variants are obtained from sequencing, the markers with structural variants can also be used in the new methods. Although the new methods have many advantages over the existing methods, they have some limitations. First, in the additive-dominant model, the numbers of marker genotypes in the two extreme pools need to be predicted when dQTG-seq1 is used, and its predicted accuracy decreases with increasing dominance ratio; in particular, when a QTG is extremely over-dominant (jd =aj > 4.0), dQTG-seq1 does not work. To address this issue, each extreme plant is sequenced, and dQTG-seq2 is used. Second, the new and existing BSA methods are affected by QTL size, sample size, sampling fraction, and sequencing depth. The best values for these factors can improve the accuracy of these methods (Zou et al., 2016) (Supplemental Table 9). In this study, we also determined the best sampling plan for the new methods ( Figure 5).
Although marker segregation distortion is a common phenomenon in bi-parental populations, only three markers around 27 known genes (Supplemental Table 6) had segregation distortion, indicating that marker segregation distortion had little importance in the analysis of real rice data. In our simulation study, marker segregation distortion affected the detection of genes (Supplemental Table 10), and its methodological extension needs to be further addressed in the near future.

Materials
The maize dataset was derived from Zhang et al. (2019), in which the F 1 plants between HZS (an elite inbred line) and 1462 (a temperate inbred line) were self-pollinated and backcrossed to HZS to produce F 2 and BC 1 F 1 seeds, respectively. The 20% tallest (580) and 20% shortest (567) plants were selected from each BC 1 F 1:2 family and bulked in equal ratios into high and low pools. The two pools were collectively subjected to whole-genome sequencing with up to >2803 coverage, and 197 021 high-quality SNPs were generated.
The genotypic and phenotypic observations from a rice IMF 2 population were downloaded from Zhou et al. (2012) (http://www.pnas.org/content/ suppl/2012/09/07/1214141109.DCSupplemental). In the IMF 2 population, there were 278 individuals and 1619 bin markers, and yield per plant in 1998 was re-analyzed. All the IMF 2 individuals were sorted according to their phenotypic values, and 20% of the extreme high and low individuals were selected to form the high and low pools. Using dQTG-seq1, the observed numbers of marker alleles were used to predict the numbers of marker genotypes. The observed and predicted numbers in extreme pools were used to calculate G w . Using dQTG-seq2, the observed numbers of marker alleles and genotypes in extreme pools were used to calculate G w .

Idea for constructing the statistic in the new method
In Figure 6A-6D, we drew the mixture distributions of one simulated QTG with four types of genetic model (a s 0 and d = 0; a = d; d = 6a; a = 0 and d s 0) in the F 2 population. We found that the difference in allele frequencies between high and low pools decreased with increasing jd =aj; in particular, no difference was observed when a = 0 and d s 0.
The theoretical derivation is shown in Supplemental File 1, and the frequency changes in marker alleles and genotypes are shown in Figure 6E. Almost all the statistics in the existing BSA methods are derived from the differences between allelic read frequencies in the two pools. This demonstrates the inability of existing BSA methods to identify QTGs with fully over-dominant effects . It should be noted that the differences in marker genotype frequencies between the two pools exist under all types of genetic model. Thus, genotypic frequencies should be used to detect all types of QTGs, including fully over-dominant QTLs. The values of G 1 , G 2, and G w are also shown in Figure 6F. As the dominance ratio increases, the value of G 1 calculated from the frequencies of marker alleles decreases, and the value of G 2 calculated from the frequencies of marker genotypes increases. In the extreme case (full over-dominance), the value of G 1 is close to zero, and the value of G 2 is close to G W . This is why the new statistic G W in this study is constructed from the frequencies of (read counts of) marker alleles and genotypes.
To estimate the five parameters above, the five equations below are used (Darvasi and Soller, 1992;Lebowitz et al., 1987).
where f ML (f QL ) and f MH (f QH ) are the frequencies of read counts of marker allele M (Q) at the currently detected marker (or QTG) in low and high pools, respectively. The five parameters above are solved using the R package BB. The initial value for jf QL À f QH j is set as the maximum difference of marker allelic read frequencies in the current chromosome between high and low pools; the initial values for x H and x L were determined based on P H , P L , a = 0, and d = 0; the initial values for a and d are set to zero; the initial value for r is set to the recombinant fraction between the putative QTL and the marker with the maximum allele frequency difference on the chromosome based on the genetic relationship between physical and genetic distances (Supplemental Table 11) and the Haldane function (Haldane, 1919). If it does not converge in the estimation of five parameters, let G 2 be zero.

Construction of a new statistic, G w
A new statistic, G w, is a weighted mean of G 1 and G 2 and is indicated by are standard G statistics (Magwene et al., 2011); E(n Qij ) and E(n Gij ) are the theoretical numbers of read counts of marker alleles and genotypes, which are calculated from the 2 3 2 and 2 3 3 contingency tables, respectively. The Smooth_G w value is a weighted average of the original statistic (G w ) values for all the SNPs in a window (Magwene et al., 2011) and is indicated by Although G 1 and G 2 obey the c 2 distribution, they are correlated. Thus, the distribution of G w is unclear. In this case, permutation experiments are used to determine the critical value of G w for significant QTGs (Fisher, 1935;Pitman, 1937).

Monte Carlo simulation studies
The purposes of a series of Monte Carlo simulation studies were to confirm whether the new methods could identify fully and extremely over-dominant and small-effect genes, to explore the conditions under which dQTG-seq1 was feasible, and to obtain the best sampling plans for the new methods.
In the second simulation experiment, one small-effect (5%) QTL was simulated under extremely (d/a = 4) and fully (a = 0 and d s 0) overdominant models to imitate a secondary F 2 population with a sample size of 500. The others were the same as those in simulated dataset I. This was simulated dataset II. In the third simulation experiment, five QTLs were simulated under extremely (jd =aj = 4) and fully (a = 0 and d s 0) over-dominant models to imitate a primary F 2 population with a sample size of 500. Their positions and effects are shown in Supplemental Table 3. The others were the same as those in simulated dataset I. This was simulated dataset III.
In the fourth simulation experiment, two QTLs with jd =aj = 0.5 were simulated on one chromosome. The sample size in the F 2 was 500, each QTL size was 10%, and QTL positions were placed at 25.00 cM and 45.00 cM. The others were the same as those in simulated dataset I. This was simulated dataset IV. In the fifth simulation experiment, we investigated the effect of sample size, QTL size, and number of plants in each extreme pool on the statistical power of dQTG-seq1 and dQTG-seq2. QTL sizes were set to 2%, 5%, and 10%; sample sizes in the F 2 were set to 250, 500, 1 000, and 2 000; and the number of plants in each extreme pool was set to 30, 50, and 80. The others were the same as those in simulated dataset I. The simulated datasets of the 36 parameter combinations were analyzed by dQTG-seq1 and dQTG-seq2 to obtain suitable sampling plans for the two new methods. This was simulated dataset V. In the sixth simulation experiment, to investigate the effect of marker segregation distortion on QTL detection, one QTL with jd =aj = 1.0 was simulated in the F 2 . The sample size was 500, the QTL position overlapped with a segregation distortion locus (SDL) at 50.00 cM, the sizes of both QTL and SDL were 10%, and the genetic models of the SDL were set to additive (a s 0 and d = 0), completely dominant (d/a = 1.0), and fully over-dominant (a = 0 and d s0). The others were the same as those in simulated dataset I. This was simulated dataset VI.
The above simulated datasets are found in Supplemental Data 1.

SUPPLEMENTAL INFORMATION
Supplemental Information is available at Plant Communications Online.