QTL Mapping and Candidate Gene Identification of Swollen Root Formation in Turnip

The swollen root is an important agronomic trait and is a determinant of yield for turnips, which are cultivated as both vegetables and fodder. However, the genetic mechanism of swollen root formation is poorly understood. In this study, we analyzed the F2 and BC1P2 populations derived from a cross between “10601” (European turnip with swollen root, Brassica rapa ssp. rapifera, AA, 2n = 2× = 20) and “10603” (Chinese cabbage with normal root, Brassica rapa ssp. pekinensis, AA, 2n = 2× = 20), and suggested that the swollen root is a quantitative trait. Two major quantitative trait loci (QTLs), FR1.1 (Fleshy root 1.1) and FR7.1 (Fleshy root 7.1), were identified by QTL-seq analysis and further confirmed by QTL mapping in F2 and BC1P2 populations. The QTL FR1.1 with a likelihood of odd (LOD) of 7.01 explained 17.2% of the total phenotypic variations for root diameter and the QTL FR7.1 explained 23.0% (LOD = 9.38) and 31.0% (LOD = 13.27) of the total phenotypic variations in root diameter and root weight, respectively. After a recombinant screening, the major QTL FR7.1 was further narrowed down to a 220 kb region containing 47 putative genes. A candidate gene, Bra003652, which is a homolog of AT1G78240 that plays an essential role in cell adhesion and disorganized tumor-like formation in Arabidopsis thaliana, was identified in this region. In addition, expression and parental allele analysis supported that Bra003652 was a possible candidate gene of QTL FR7.1 for swollen root formation in turnip. Our research may provide new insight into the molecular mechanism of swollen root formation in root crops.


Introduction
Turnip (Brassica rapa ssp. rapifera, AA, 2n = 2× = 20), which is cultivated for use as a vegetable, medicine, and fodder, is characterized by a swollen root. The turnip represents an important type in B. rapa, which displays considerable morphological diversity (leafy vegetables, oil crop, and stem/root crop). The swollen root is an important storage organ, which contains sugar, proteins, minerals, and many kinds of vitamins, and determines the yield and quality of turnip. Therefore, understanding the molecular mechanism underlying swollen root formation is important for improving root crops performance. Some researchers have shown that the swollen root of turnip is a taproot [1][2][3], while a few studies have mentioned that the thickened part of turnips consists of both hypocotyl and root [4][5][6]. Anatomical observation has identified that morphological changes of the xylem lead to an initiation of swollen root formation [6]. However, the genetic mechanism underlying swollen root formation and its regulation remains unknown. With respect to turnips, the root shape and size are quantitative traits [7]. A previous study has revealed that a major quantitative trait loci (QTL), TuQTL-1 near BrFLC2 (Flowering Locus C), and two novel QTLs on chromosomes A1 and A5 are involved in swollen root formation [8]. In addition, 18 QTLs controlling the traits of fleshy roots, including seven QTLs for root width, five QTLs for root length, and six QTLs for root weight, have been identified based on an F 2 population [3]. However, no swollen root-related QTLs have been fine-mapped thus far, and this makes it difficult to perform for further research.
The de novo assembly of the Chinese cabbage (Chiifu-401-42) [9] genome greatly promotes the development of functional genomics of Chinese cabbage and other Brassica crops. More and more genomic and transcriptomic studies have been used to study the swollen root traits [10]. The analysis of 199 B. rapa and 119 B. oleracea accessions has identified a large number of genes related to the synthesis, transport, or binding of sugar or cellulose, the expansion gene family, root development, and cell division, which participates in the formation of the storage root in turnip [11]. High-throughput analysis of small ribonucleic acid (RNA) sequences has indicated that various microRNAs may play important roles in the growth, development, and response to dark environments of turnip [12]. Comparative transcriptomics has revealed 841 differentially expressed genes (DEGs), which are homologues in turnips and kohlrabi [13]. The Bra-FLOR1 and Bra-CYP735A2 genes are involved in tuber initiation and/or tuber growth based on cytological, physiological, genetic, and transcriptomic approaches [6]. However, the function of these genes and microRNAs needs further verification.
The genetics underlying storage organ formation in many other crops can be used to provide a reference for turnip analysis. Species with storage organs, including radish, potato, rutabaga, carrot, beets, sweet potato, and cassava, play significant roles in agriculture, and their storage organs are controlled by complex interactions involving genetic, environmental, and physiological factors [5,6,14]. Sucrose is crucial in the development of storage organs in many root/stem crops [15,16]. The sucrose metabolism pathways are the most significantly activated in radish storage roots, and the expression levels of sucrose synthase 1 (SUS1) are correlated with root thickening rates [14,17], which is consistent with the study that sucrose synthase is a major enzyme involved in the storage root development of radish [16]. In addition, plant hormones, such as cytokinin, auxin, jasmonic acid (JA), abscisic acid (ABA), and gibberellin (GA) are important regulatory factors in the formation and thickening of tuberous roots [18][19][20][21][22]. Furthermore, short days and cool temperatures promote tuber formation in potato [23][24][25].
In this study, we observed the development and inheritance pattern of the swollen root traits in the progeny of a cross between two subspecies of B. rapa, European turnip with swollen root and Chinese cabbage with normal root. QTL-seq technology was used to sequence the deoxyribonucleic acid (DNA) samples from large root and small root pools and to rapidly identify major QTLs associated with swollen root traits. Then, two major QTLs were confirmed and fine-mapped using classical linkage mapping. Expression pattern analyses and comparative genomics were conducted for the candidate genes. Our research provides some candidate genes and lays a foundation for studying the molecular mechanism of swollen root formation.

Inheritance of Swollen Root in B. rapa
Significant differences in root-related traits were observed among the "10601" (P1, swollen root), "10603" (P2, normal root), and F 1 hybrids, in 2018 and 2019 ( Figure 1A and Table 1). P1 exhibited a larger swollen root than P2, with an average root diameter of 6.4-7.0 cm, root weight of 316.0-330.0 g, and root length of 13.5-14.1 cm; P2 had a small root with an average root diameter of 1.5-1.9 cm, root weight of 23.8-30.8 g, and root length of 2.5-12.0 cm, in 2018 and 2019 (Table 1). The value for the F 1 hybrid (except root length in 2019) was intermediate between those of its two parents, but much closer to that of P1 (Table 1). F 1 plants derived from reciprocal crosses between P1 and P2 exhibited obvious swollen roots, which suggested that the swollen root trait was dominant in this cross ( Figure 1A). length in 2019) was intermediate between those of its two parents, but much closer to that of P1 (Table 1). F1 plants derived from reciprocal crosses between P1 and P2 exhibited obvious swollen roots, which suggested that the swollen root trait was dominant in this cross ( Figure 1A).
The root diameter, root weight, and root length traits showed large coefficients of variation (CV%) among the F2 and BC1P2 populations, and continuous distribution in three segregating populations (F2 in 2018, F2 in 2019, and BC1P2 in 2019), which indicated that the root-related traits are quantitative in B. rapa ( Figure 1B and Table 1).  The root diameter, root weight, and root length traits showed large coefficients of variation (CV%) among the F 2 and BC 1 P 2 populations, and continuous distribution in three segregating populations (F 2 in 2018, F 2 in 2019, and BC 1 P 2 in 2019), which indicated that the root-related traits are quantitative in B. rapa ( Figure 1B and Table 1).

Phenotypic and Anatomical Observation of the Hypocotyl and Root from Parents
Phenotypic observation of parent plants in different development stages showed that the swollen root of turnip consisted of hypocotyl and root tissues ( Figure 2). The diameter of the hypocotyl showed an increasing trend in P1 and P2 from 21 to 49 days after germination (DAG) (Figure 2A,D). The root had an obvious cortex splitting phenotype on 35 DAG (Figure 2A, red arrow) and the root diameter exhibited a significant increase from 35 DAG in P1 plants, which indicated that the cortex splitting stage was also a sign of a transition from primary growth to secondary growth in turnip, similar to radish [17,26,27].

Phenotypic and Anatomical Observation of the Hypocotyl and Root from Parents
Phenotypic observation of parent plants in different development stages showed that the swollen root of turnip consisted of hypocotyl and root tissues ( Figure 2). The diameter of the hypocotyl showed an increasing trend in P1 and P2 from 21 to 49 days after germination (DAG) (Figure 2A,D). The root had an obvious cortex splitting phenotype on 35 DAG (Figure 2A, red arrow) and the root diameter exhibited a significant increase from 35 DAG in P1 plants, which indicated that the cortex splitting stage was also a sign of a transition from primary growth to secondary growth in turnip, similar to radish [17,26,27].
The bottom part of hypocotyls from the two parents was selected for the anatomical sections, and we found that the phloem (Ph), vascular cambium (VC), xylem (X), vessel (VE), and pith (P) could be clearly distinguished among the five developmental stages ( Figure 2B). From 21 to 35 DAG, the xylem width of P1 was narrower than that of P2. Then, the xylem width increased significantly from 35 to 56 DAG in P1, while that of P2 showed limited radial growth ( Figure 2B,C,E,F). There was little difference in the phloem between P1 and P2. These results indicated that the hypocotyl of turnip strongly increased in diameter, mainly by the rapidly expanding xylem parenchyma layer. The bottom part of hypocotyls from the two parents was selected for the anatomical sections, and we found that the phloem (Ph), vascular cambium (VC), xylem (X), vessel (VE), and pith (P) could be clearly distinguished among the five developmental stages ( Figure 2B). From 21 to 35 DAG, the xylem width of P1 was narrower than that of P2. Then, the xylem width increased significantly from 35 to 56 DAG in P1, while that of P2 showed limited radial growth ( Figure 2B,C,E,F). There was little difference in the phloem between P1 and P2. These results indicated that the hypocotyl of turnip strongly increased in diameter, mainly by the rapidly expanding xylem parenchyma layer.

Quantitative Trait Loci (QTL)-seq Predicted Candidate Regions for Controlling the Swollen Root Traits
In total, 15.47-15.99 Gb clean reads were obtained for the P1, P2, L-pool, and S-pool. The average sequencing depth ranged from 29.21 to 33.86×, and the coverage ranged from 89.79 to 95.76% (Table S1). Using the B. rapa genome V1.5 as reference, 100,565-261,027 insertion and deletion (InDel) and 427,823-1,369,893 single nucleotide polymorphisms (SNP) variations were identified for P1, P2, L-pool, and S-pool (Table S2). Finally, 753,153 high-quality SNPs that were homozygous for each parent and polymorphic in the two parents were obtained and used for calculating the SNP-index and ∆ (SNP-index). A distribution map of the SNPs was drawn by plotting the number of SNPs using a sliding window of 100 kb ( Figure 3).

Quantitative Trait Loci (QTL)-seq Predicted Candidate Regions for Controlling the Swollen Root Traits
In total, 15.47-15.99 Gb clean reads were obtained for the P1, P2, L-pool, and S-pool. The average sequencing depth ranged from 29.21 to 33.86×, and the coverage ranged from 89.79 to 95.76% (Table S1). Using the B. rapa genome V1.5 as reference, 100,565-261,027 insertion and deletion (InDel) and 427,823-1,369,893 single nucleotide polymorphisms (SNP) variations were identified for P1, P2, L-pool, and S-pool (Table  S2). Finally, 753,153 high-quality SNPs that were homozygous for each parent and polymorphic in the two parents were obtained and used for calculating the SNP-index and ∆ (SNP-index). A distribution map of the SNPs was drawn by plotting the number of SNPs using a sliding window of 100 kb ( Figure 3).

Confirmation and Fine-Mapping of the Two QTLs: FR1.1 and FR7.1
Traditional QTL analysis, based on 165 F2 plants and 153 BC1P2 plants in two environments, was performed to verify the QTLs FR1.1 and FR7.1 detected by QTL-seq analysis. Fifty-two pairs of InDel markers and 39 pairs of SNP markers at the target regions of the two QTLs were designed, based on the whole genome sequencing data of the parents (Tables S3 and S4).
Among the developed markers on chromosome 7, 13 markers that covered the QTL region and showed polymorphism between parents were used to genotype the segregating population for QTL analysis. In the BC1P2 population, FR7.1 was narrowed to the 17.01-17.80 Mb region between the SNP markers A07_S04 and A07_S21, which explained 13.10% and 16.00% of the phenotypic variation of root diameter and root weight, respectively (Table 2 and Figure 5A). In the F2 population, FR7.1 was located within the 17.31-17.60 Mb region between the SNP markers A07_S43 and A07_S06 and explained 23.00% and 31.00% of the phenotypic variation of root diameter and root weight with likelihood of odd (LOD) values of 9.38 and 13.27, respectively (Table 2 and Figure 5A). By combining the QTL mapping results from the F2 and BC1P2 populations, the QTL FR7.1 was confirmed and mapped to the marker interval between A07_S43 (17.31 Mb) and A07_S06 (17.60 Mb) on chromosome 7. After a recombinant screening using 716 plants from F2-2018, F2-2019, and BC1P2-2019 populations, the QTL FR7.1 was further narrowed down to a 220 kb region (17.31-17.53 Mb) bounded by A07_S43 and A07_S35 ( Figure  5B,C).   Traditional QTL analysis, based on 165 F 2 plants and 153 BC 1 P 2 plants in two environments, was performed to verify the QTLs FR1.1 and FR7.1 detected by QTL-seq analysis. Fifty-two pairs of InDel markers and 39 pairs of SNP markers at the target regions of the two QTLs were designed, based on the whole genome sequencing data of the parents (Tables S3 and S4).
Among the developed markers on chromosome 7, 13 markers that covered the QTL region and showed polymorphism between parents were used to genotype the segregating population for QTL analysis. In the BC 1 P 2 population, FR7.1 was narrowed to the 17.01-17.80 Mb region between the SNP markers A07_S04 and A07_S21, which explained 13.10% and 16.00% of the phenotypic variation of root diameter and root weight, respectively (Table 2 and Figure 5A). In the F 2 population, FR7.1 was located within the 17.31-17.60 Mb region between the SNP markers A07_S43 and A07_S06 and explained 23.00% and 31.00% of the phenotypic variation of root diameter and root weight with likelihood of odd (LOD) values of 9.38 and 13.27, respectively (Table 2 and Figure 5A). By combining the QTL mapping results from the F 2 and BC 1 P 2 populations, the QTL FR7.1 was confirmed and mapped to the marker interval between A07_S43 (17.31 Mb) and A07_S06 (17.60 Mb) on chromosome 7. After a recombinant screening using 716 plants from F 2 -2018, F 2 -2019, and BC 1 P 2 -2019 populations, the QTL FR7.1 was further narrowed down to a 220 kb region (17.31-17.53 Mb) bounded by A07_S43 and A07_S35 ( Figure 5B,C). which was physically located in the region of 4.62-5.29 Mb on chromosome 1. This QTL explained 17.2% of the phenotypic variation in root diameter with a LOD of 7.01 (Table  3). Thus, FR1.1 was also confirmed and narrowed down to a 670 kb region (4.62-5.29 Mb) on chromosome 1.  Similarly, 10 polymorphic InDel markers that covered the FR1.1 region were used for linkage analysis in the F 2 population. A major QTL for root diameter between two InDel markers, A01_ID03 and A01_ID31, was identified by QTL mapping analysis, which was physically located in the region of 4.62-5.29 Mb on chromosome 1. This QTL explained 17.2% of the phenotypic variation in root diameter with a LOD of 7.01 (Table 3). Thus, FR1.1 was also confirmed and narrowed down to a 670 kb region (4.62-5.29 Mb) on chromosome 1.

Identification of Candidate Genes for FR7.1 in Turnip
Based on the results of QTL-seq and QTL mapping, FR7.1 was selected for further analysis. The candidate interval contained 47 predicted genes in the B. rapa reference genome (http://brassicadb.org/brad/index.php). Among these genes, 31 genes have SNP variation between the two parents (Table S5). According to the current research, sucrose metabolism, hormonal and the cell cycle pathways are the central components of tuberization in radish, potato, and B. juncea [6,17,28,29]. Therefore, four genes involved in the cell cycle, glycolysis/gluconeogenesis, and plant hormone signal transduction were selected, namely, Bra003606, Bra003649, Bra003650, and Bra003652, according to the annotations ( Figure 5C and Table S6).
In order to determine whether the expression levels of the candidate genes may contribute to the variance of swollen root traits, we studied the expression patterns of Bra003606, Bra003649, Bra003650, and Bra003652 in the hypocotyls of P1 and P2 via quantitative real-time polymerase chain reaction (qRT-PCR). The results showed that the expression level of Bra003652 was significantly higher in P1 than in P2 (p < 0.01). However, the expression levels of Bra003606, Bra003649, and Bra003650 did not differ between P1 and P2 ( Figure 6). Therefore, Bra003652 aroused our attention.

Identification of Candidate Genes for FR7.1 in Turnip
Based on the results of QTL-seq and QTL mapping, FR7.1 was selected for further analysis. The candidate interval contained 47 predicted genes in the B. rapa reference genome (http://brassicadb.org/brad/index.php). Among these genes, 31 genes have SNP variation between the two parents (Table S5). According to the current research, sucrose metabolism, hormonal and the cell cycle pathways are the central components of tuberization in radish, potato, and B. juncea [6,17,28,29]. Therefore, four genes involved in the cell cycle, glycolysis/gluconeogenesis, and plant hormone signal transduction were selected, namely, Bra003606, Bra003649, Bra003650, and Bra003652, according to the annotations ( Figure 5C and Table S6).
In order to determine whether the expression levels of the candidate genes may contribute to the variance of swollen root traits, we studied the expression patterns of Bra003606, Bra003649, Bra003650, and Bra003652 in the hypocotyls of P1 and P2 via quantitative real-time polymerase chain reaction (qRT-PCR). The results showed that the expression level of Bra003652 was significantly higher in P1 than in P2 (p < 0.01). However, the expression levels of Bra003606, Bra003649, and Bra003650 did not differ between P1 and P2 ( Figure 6). Therefore, Bra003652 aroused our attention. The physical location of the Bra003652 gene is 17,537,326-17,540,584 bp, which contains the SNP marker A07_S35 at which the LOD value curve peaked for FR7.1 ( Figure  5C, Tables 2 and S6). Bra003652 is a homolog of AT1G78240, which encodes the tumorous shoot development 2 (TSD2) protein and is related to root development in Arabidopsis thaliana [30]. Bra003652 and AT1G78240 are 85.38% identical at the nucleotide level (Figure S1) and share 86.14% amino acid identity ( Figure S2). This suggests that these genes encode enzymes with similar functions.
To define the structure of Bra003652, the amplified products of the complementary DNA (cDNA) and genomic DNA (gDNA) fragments were sequenced using Bra003652-specific primers. The result showed that Bra003652 contains eight exons and seven introns, encoding a protein of 684 amino acids (Figures 7A and S2). The putative protein contains a putative S-adenosyl-L-methionine-dependent methyltransferase do- The physical location of the Bra003652 gene is 17,537,326-17,540,584 bp, which contains the SNP marker A07_S35 at which the LOD value curve peaked for FR7.1 ( Figure 5C, Table 2 and Table S6). Bra003652 is a homolog of AT1G78240, which encodes the tumorous shoot development 2 (TSD2) protein and is related to root development in Arabidopsis thaliana [30]. Bra003652 and AT1G78240 are 85.38% identical at the nucleotide level ( Figure S1) and share 86.14% amino acid identity ( Figure S2). This suggests that these genes encode enzymes with similar functions.
To define the structure of Bra003652, the amplified products of the complementary DNA (cDNA) and genomic DNA (gDNA) fragments were sequenced using Bra003652specific primers. The result showed that Bra003652 contains eight exons and seven introns, encoding a protein of 684 amino acids ( Figure 7A and Figure S2). The putative protein contains a putative S-adenosyl-L-methionine-dependent methyltransferase domain ( Figure S2).
Comparative sequence analysis was used to align the genomic DNA sequences of Bra003652 from both parents. Some variations were detected in the open reading frame of Bra003652, including thirteen SNP variations and a three-base deletion between the two parents ( Figure 7A). A multiple protein sequence alignment of Bra003652 and its homologous genes from five inbred lines, including two materials with swollen root (P1 and radish) and three materials with normal root (Arabidopsis, Chiifu, and P2), was conducted. There were three amino acid variations and one amino acid deletion between P1 and P2, which did not occur in the region of the putative domain ( Figure 7B and Figure S2). Among the four variations, the first three variations (S → N, an amino acid deletion and S → P) were the same for P1 and radish, and they were different among Arabidopsis, Chiifu, and P2 ( Figure 7B, black arrows). This gave rise to a hypothesis that these variations might influence the function of Bra003652 and involve in swollen root formation in turnips.  Figure S2). Comparative sequence analysis was used to align the genomic DNA sequences of Bra003652 from both parents. Some variations were detected in the open reading frame of Bra003652, including thirteen SNP variations and a three-base deletion between the two parents ( Figure 7A). A multiple protein sequence alignment of Bra003652 and its homologous genes from five inbred lines, including two materials with swollen root (P1 and radish) and three materials with normal root (Arabidopsis, Chiifu, and P2), was conducted. There were three amino acid variations and one amino acid deletion between P1 and P2, which did not occur in the region of the putative domain ( Figures 7B and S2). Among the four variations, the first three variations (S → N, an amino acid deletion and S → P) were the same for P1 and radish, and they were different among Arabidopsis, Chiifu, and P2 ( Figure 7B, black arrows). This gave rise to a hypothesis that these variations might influence the function of Bra003652 and involve in swollen root formation in turnips. At the same time, we also analyzed the sequences of Bra003606, Bra003649, and Bra003650 from parents. Some variations were found in the coding sequences of Bra003649 and Bra003606 (Figures S3 and S5), while that of parents were identical in Bra003650. The putative protein sequences of parents were identical in Bra003649 ( Figure  S4), while an amino acid variation was found between P1 and P2 in Bra003606 ( Figure  S6). Further comparative analysis showed that the mutation type of this amino acid was consistent in P1 (turnip) and Chiifu (Chinese cabbage), but different from P2 (Chinese cabbage), which indicated that the variation of Bra003606 might not be involved in the formation of swollen root ( Figure S6). These results further support that Bra003652 is the most possible candidate gene in FR7.1

QTL-seq Technology Accelerated the Research of Root-Related Traits
Genetic analysis shows that most root-related traits are quantitative traits controlled by multiple genes [3,7]. In our study, a European turnip and a Chinese cabbage were selected as the parental lines and used to construct the F1, F2, and BC1P2 populations. Genetic analysis showed that these traits exhibited quantitative inheritance patterns, consistent with the previous research results (Figure 1).
Many QTLs related to root traits have been identified and are scattered on multiple chromosomes in turnip, but a few candidate genes have been identified [3,7,8]. Traditional QTL mapping includes primary mapping, fine mapping, and the identification of candidate genes, which are time-consuming and expensive. QTL-seq technology, which has been applied to the QTL mapping of many traits in different species [31][32][33], was used and identified two QTLs, FR7.1 and FR1.1, in this research (Figure 4). Subsequently, traditional QTL analysis of the F2 and BC1P2 populations was used to verify the accuracy of QTL-seq analysis. InDel and SNP markers were developed and used to construct linkage maps in segregated populations. InDel is widely used in traditional QTL map- At the same time, we also analyzed the sequences of Bra003606, Bra003649, and Bra003650 from parents. Some variations were found in the coding sequences of Bra003649 and Bra003606 (Figures S3 and S5), while that of parents were identical in Bra003650. The putative protein sequences of parents were identical in Bra003649 ( Figure S4), while an amino acid variation was found between P1 and P2 in Bra003606 ( Figure S6). Further comparative analysis showed that the mutation type of this amino acid was consistent in P1 (turnip) and Chiifu (Chinese cabbage), but different from P2 (Chinese cabbage), which indicated that the variation of Bra003606 might not be involved in the formation of swollen root ( Figure S6). These results further support that Bra003652 is the most possible candidate gene in FR7.1.

QTL-seq Technology Accelerated the Research of Root-Related Traits
Genetic analysis shows that most root-related traits are quantitative traits controlled by multiple genes [3,7]. In our study, a European turnip and a Chinese cabbage were selected as the parental lines and used to construct the F 1 , F 2 , and BC 1 P 2 populations. Genetic analysis showed that these traits exhibited quantitative inheritance patterns, consistent with the previous research results (Figure 1).
Many QTLs related to root traits have been identified and are scattered on multiple chromosomes in turnip, but a few candidate genes have been identified [3,7,8]. Traditional QTL mapping includes primary mapping, fine mapping, and the identification of candidate genes, which are time-consuming and expensive. QTL-seq technology, which has been applied to the QTL mapping of many traits in different species [31][32][33], was used and identified two QTLs, FR7.1 and FR1.1, in this research ( Figure 4). Subsequently, traditional QTL analysis of the F 2 and BC 1 P 2 populations was used to verify the accuracy of QTL-seq analysis. InDel and SNP markers were developed and used to construct linkage maps in segregated populations. InDel is widely used in traditional QTL mapping by polyacrylamide gel electrophoresis (PAGE) or agarose gel electrophoresis technology, which is accurate but time-consuming and labor-intensive. With the development of next-generation sequencing technologies, the genome of B. rapa was sequenced and assembled in 2011 [9]. Subsequently, many materials with different ecological types have been sequenced or re-sequenced by whole genome sequencing (WGS) or specific length-amplified fragment sequencing [11,34], which has revealed a large number of SNPs among different materials. In our research, we obtained 753,153 high-quality SNPs that showed polymorphism between P1 and P2. In addition, SNPs can be genotyped by Kompetitive Allele Specific PCR (KASP), which is one of the SNP genotyping platforms and has evolved to be global benchmarking technology [35]. SNP markers were chosen to construct a local linkage map and to fine-map FR7.1. Then, FR7.1 was narrowed to a 220 kb region between two SNP markers, A07_S04 and A07_S35 (Table 2 and Figure 5), consistent with the result of QTL-seq analysis. The results showed that the method of combining QTL-seq and QTL mapping can rapidly and accurately identify QTL.

The Authenticity and Reliability of FR1.1 and FR7.1
Previously, Cheng has shown that the domestication of tuber-forming morphotypes was the result of parallel selection and convergent crop domestication and identified 24 and 26 genomic regions under selection in turnip using reduction of diversity (ROD) and PiHS population-based integrated haplotype score (PiHS), respectively [11]. Some of these identified genomic regions have overlap with the target regions of FR1.1 and FR7.1, which further proves the reliability of FR1.1 and FR7.1. However, the two QTLs have not been reported by other researchers through traditional genetics and may be the new loci involved in the formation of swollen root.
In our study, FR7.1 was narrowed to 220 kb between two SNP markers, A07_S45 and A07_S35, which explained 23% and 31% of the phenotypic variation in root diameter and root weight in the F 2 population in 2018, respectively (Table 2). Then, we used the F 2 population and BC 1 P 2 population to verify the authenticity of FR7.1, in 2019. However, we were unable to scan any QTL in the 16.01-21.99 Mb on chromosome 7 in the F 2 population (data not shown). In 2019, many plants from F 2 populations were infected with clubroot or soft rot diseases, which affected the measurements of root-related traits. Fortunately, FR7.1 was localized to a 790 kb physical interval by the two markers A07_S04 and A07_S21, contributing 13.1% and 16% of the phenotypic variation in root diameter and root weight, respectively, in the BC 1 P 2 populations, which overlapped with the result of QTL mapping from 2018 (Table 2 and Figure 5). The result showed the QTL FR7.1 related to root diameter and root weight is reliable and could be detected in 2018 and 2019.
We also narrowed the region of the FR1.1 locus to 670 kb between A01_ID03 and A01_ID31. FR1.1 explained 17.2% of phenotypic variation in root diameter in the F 2 population, in 2018 (Table 3). The locus needed to be further verified and fine-mapped using new segregating populations, which might provide new genetic resources for the research of swollen root formation.

Analysis and Determination of Candidate Genes
In total, 47 genes were identified in the candidate region of FR7.1. On the basis of the Previous studies shows that plant hormones and sucrose are important regulators of root formation and development [2,16,18,19,36]. In addition, environmental factors including short sunshine and low temperatures can also promote tuber formation in potato [23][24][25]. Among 47 genes, four genes involved in the cell cycle, glycolysis/gluconeogenesis, plant hormone signal transduction, and root development were identified, namely Bra003606, Bra003649, Bra003650, and Bra003652 (Table S6). According to the anatomical observation, the number and size of xylem cells influenced the width of the xylem and are related to the enlargement of the hypocotyl in turnip (Figure 2). The expression models for Bra003606, Bra003649, Bra003650, and Bra003652 that were obtained via qRT-PCR support that Bra003652 is the most likely candidate gene of FR7.1 in B. rapa ( Figure 6).
Bra003652 is a homolog of AT1G78240, which encodes the TSD2 protein in A. thaliana. Mutation of the TSD2 gene reduces cell adhesion and strongly affects individuals, resulting in non-coordinated shoot development, which leads to disorganized tumor-like growth in A. thaliana [30]. Some studies have shown that the TSD gene depended on cytokinin to negatively regulate the activity of the meristem during the development of Arabidopsis [34]. In rice, the OsTSD2 regulates the root development and cellular adhesion by affecting pectin synthesis [37]. Cheng [11] has obtained a large number of genes related to the formation of swollen root in turnip, including Bra003652.
A multiple sequence alignment analysis identified three amino acid variations (S → N, S → P and S → N) and an amino acids deletion between P1 and P2 ( Figure 7B). Among these four variations, three amino acid variations were the same in the protein sequences of Bra003652 and Rsa10029866 from P1 and radish (swollen root), and these were different from the variations observed in Chiifu, P2, and Arabidopsis (normal root, Figure 7B, black arrows). These results gave rise to a hypothesis that these variations might influence the function of Bra003652 and be involved in the swollen root formation in turnips.

Root-Related Traits Positively Affect the Yield of Root Crops
Swollen root is the decisive factor for the yield of root-type crops and consists of the following three important indexes: root diameter, root length, and root weight. In addition, many above-ground traits, including plant height, plant width, and leaf weight, directly or indirectly impact the yield of root/stem crops. Correlation analysis has shown that root yield of radish, sweet potato, and carrot was positively correlated with root length and root diameter [38,39]. In our research, the correlation coefficient for root weight and root diameter was 0.71, which was higher than other correlation coefficients in the F 2 population ( Figure S7). The correlation coefficient between root weight and root length was 0.55, while the correlation coefficient between root length and root diameter was 0.47 ( Figure S7). The result shows that root diameter and root length are positively correlated with root weight, which is consistent with previous conclusions. In addition, root weight is related to plant height, plant width, and leaf weight; the correlation coefficients are 0.48, 0.54 and 0.51, respectively ( Figure S7). Root diameter has a greater effect on root weight than other traits from the above result. In the study, FR1.1 explained 17.2% phenotypic variation in root diameter and FR7.1 explained 23% phenotypic variation of root diameter, which may harbor candidate genes involved in root diameter development. Our research may provide some candidate genes and lay a foundation for improving the yield of turnip.

Parent Materials and Population Construction
Two inbred lines, Europe turnip "10601" (P1, swollen root) and Chinese Cabbage "10603" (P2, normal root) from B. rapa, both were collected by the Institute of Vegetables and Flowers, Chinese Academy of Agricultural Sciences, and were used as parents in a reciprocal cross (P1 × P2 and P2 × P1) to produce F 1 lines. F 1 (P1 × P2) was self-pollinated and backcrossed to generate the F 2 and BC 1 P 2 populations. Finally, the following five populations were obtained: P1, P2, F 1 (P1 × P2), F 1 (P2 × P1), and F 2 and BC 1 P 2 .

Phenotype Assessment
The parental line, and the F 1 and F 2 populations were seeded in commercial potting mix in 54 × 28 cm flats with 32-wells in a solar greenhouse on 27 July 2018. Then, the seedlings were transplanted to an experimental farm on 16 August 2018. The row spacing and plant spacing were 50 and 40 cm, respectively. On 6 August 2019, the parental lines, and the F 1 , F 2 , and BC 1 P 2 populations were seeded directly at the experimental farm. The row spacing and plant spacing were 60 and 50 cm, respectively. In 2018 and 2019, all experiments were carried out under standard field conditions in the Institute of Vegetables and Flowers, Chinese Academy of Agricultural Sciences (Shunyi, Beijing, China), with watering and fertilizing occurring in a timely manner. Three types of fertilizer, including organic fertilizer, compound fertilizer (N/P/K = 15:15:15) (Luxi Chemical Group Co., Ltd., Liaocheng, China), and punching fertilizer (N/P/K = 30:9:12) (Linong Feng Biotechnology Co., Ltd., Beijing, China), were used in the experimental field. The organic fertilizer including cattle manure (10 m 3 /667 m 2 ) and chicken manure (2 m 3 /667 m 2 ) and compound fertilizer (50 kg/667 m 2 ) were applied before sowing. The punching fertilizer (5-10 kg/667 m 2 ) was applied five times during the growth and development of the materials.
There were 15, 15, 15, and 447 plants in P1, P2, F 1 , and F 2 populations in 2018, respectively. In 2019, there were 15, 15, 15, 278, and 275 plants in P1, P2, F 1 , F 2 and BC 1 P 2 , respectively. For each plant growing 90 days after germination, six traits including root diameter, root length, root weight, plant height, plant width, and leaf weight were measured and recorded with a ruler and electronic scale. The root weight was measured after cleaning and drying. The root diameter and root length data were recorded by measuring the transverse diameters and largest longitudinal length of the roots, respectively. The correlation coefficients among different traits were calculated by Pearson correlation method and drawn with the R package "corrplot". The statistical analyses were performed using the SAS program.

QTL-seq Analysis
According to the phenotypic investigation of the F 2 population in the winter of 2018, 30 individuals with extremely large roots and 30 individuals with extremely small roots were selected to build a large-root DNA pool (L-pool) and a small-root DNA pool (S-pool) by mixing the same amounts of DNA (Table S1). The genomic DNAs of "10601" and "10603" were obtained from a single plant. Genomic DNA was isolated from fresh leaves using a modified cetyltrimethyl ammonium bromide (CTAB) method [40]. The two parents and two pools were sequenced using an Illumina HiSeq 2500 instrument (San Diego, CA, USA) with 150 bp paired-end sequencing reads. The average sequence depth exceeded 30× genome coverage for each sample. The adaptor reads, unknown sequences "N" (reads containing unknown nucleotides >10%) and low-quality reads (reads containing more than 50% bases with a Q value ≤ 3) were removed from the raw data using the NGSQC toolkit [41]. The clean reads from the two extreme pools and two parents were mapped to the B. rapa genome V1.5 (http://brassicadb.org/brad/datasets/pub/BrassicaceaeGenome/ Brassica_rapa/V1.0/Bra_Chromosome_V1.5/) using Burrows-Wheeler Aligner (BWA) software with the BWA-MEM algorithm [42]. Based on the sequence alignment results, SNPs and InDels were called using Samtools and GATK (The Genome Analysis Toolkit) software. The density distribution map of the SNPs was drawn with a sliding window with a physical distance of 100 kb. Then, Annovar software was used to annotate the SNPs and InDels. Following published methods [33,43], the SNP-index for each SNP position was calculated.
The potential candidate region (p < 0.05) was determined by the ∆ (SNP-index) curve. P2 was selected as a reference to calculate the SNP-index and ∆ (SNP-index) of the two bulk pools according to the Takagi method [33]. An SNP-index = 0.5 means that both parents have an equal contribution to the bulked progeny. If an SNP locus is related to the target trait, the SNP-index of that locus will deviate significantly from 0.5. An SNP-index = 0 indicates that all short reads are from P2. Conversely, an SNP-index = 1 indicates that all short reads are from P1. Then, we calculated the average value in a 1 Mb sliding window based on the SNP-index and ∆ (SNP-index) values of each SNP locus, and drew a graph showing the 95% confidence interval. If the SNPs in the two extreme pools are consistent with other traits except for the primary target region, the regions related to the target trait will show linkage disequilibrium. The area above the threshold is the area where the candidate QTL is located in.

Molecular Marker Development, Linkage Analysis, and QTL Mapping
To validate the QTL-seq analysis results, traditional QTL mapping was carried out. InDel markers were designed using Primer 5.0 software, and SNP markers were designed by the Laboratory of the Government Chemist based on the B. rapa genome v1.5. All primer pairs were synthesized by Sangon Biological and Engineering Co (Shanghai, China). Linkage analysis was performed by applying the JoinMap v4.0 [44] to the Kosambi function, and the maximum-likelihood method was used to calculate the genetic distance. MapQTL v4.0 [45] was applied to detect the QTLs using the Interval Mapping (IM) and Multiple QTL Mapping (MQM) pattern under a threshold LOD = 2.0.

RNA Extraction and Quantitative Real-Time PCR
In order to analyze the candidate genes related to the swollen root traits, expression patterns of the candidate genes (Bra003606, Bra003649, Bra003650, and Bra003652) were investigated by qRT-PCR. In the spring of 2020, P1 and P2 materials were seeded in the greenhouse (Haidian, Bejing, China), with temperatures ranging from 14 • C at night to 25 • C in the daytime. Six hypocotyl samples of P1 and P2 were collected at 21,35,42,49, and 56 days after germination. Three samples were used for Paraffin section and the other three plants were used to harvest RNA. Total RNA was extracted using the Trizol reagent (Thermo Fisher Scientific, Shanghai, China), according to the manufacturer's instructions. Reverse transcription was conducted by EasyScript ® One-Step gDNA Removal and cDNA Synthesis SuperMix (TransGen Biotech, Beijing, China), using a standard protocol. Real-time PCR was performed using ChamQ Universal SYBR qPCR Master Mix (Vazyme Biotech Co., Ltd., Nanjing, China) on an ABI QuantStudio 12K Flex real-time PCR system (Applied Biosystem, Foster City, CA, USA), according to the manufacturer's instructions. The Actin gene of Chinese cabbage was used as an internal control (forward, 5 -GGAGCTGAGAGATTCCGTTG-3 , and reverse, 5 -GAACCACCACTGAGGACGAT-3 ) [46]. The special primers were designed using the Primer-BLAST tool (https://www. ncbi.nlm.nih.gov/tools/primer-blast/index.cgi?LINK_LOC = BlastHome). The primer sequences are listed in Table S7. The 2 −∆∆Ct method [47] was used to determine the relative expression levels of the genes. Then, the average relative expression levels were calculated, and t-tests were performed using the SAS program to test the significance of the differences in expression levels among the different samples.

Histological Analysis of Hypocotyls
The hypocotyls of "10601" and "10603" plants were collected at 21, 35, 42, 49, and 54 days after germination and were fixed, embedded, cut into slices, and stained according to the conventional paraffin-sectioning methods. The slices were observed and photographed using the microscope (Nikon ECLIPSE 80i, Tokyo, Japan). The xylem and phloem width were measured with Image J software (https://imagej.en.softonic.com/).

Sequence Analysis of Candidate Genes
In order to clarify the difference of candidate genes between parents, the code sequences (CDS) and gene sequence were amplified from the cDNA and gDNA of P1 and P2, respectively. The open reading frame and protein sequences of Bra003652 and Bra003606 from Chiifu were collected, according the B. rapa v1.5 genome (http://brassicadb.org/brad/ index.php). The protein sequence of AT1G78240 and Rsa10029866 which are homologous genes of Bra003652 from A. thaliana and radish were obtained, according to the corresponding reference genome [48,49]. Nucleotide and protein sequences were aligned using SnapGene software, MultAlin (http://multalin.toulouse.inra.fr/multalin/multalin.html) and MUSCLE (https://www.ebi.ac.uk/Tools/msa/muscle/).

Conclusions
In summary, anatomical observation indicated that the swollen root of turnip consists of both hypocotyl and root tissues and the swollen root traits were quantitative. Two major QTLs, FR1.1 and FR7.1 involved in the root diameter and root weight, were detected by QTL-seq technology and classical QTL mapping in the turnip × Chinese cabbage population. The major QTL, FR7.1, explained 31% (LOD = 13.27) and 23.00% (LOD = 9.38) of the total phenotypic variations in root weight and root diameter, respectively. Expression pattern and comparative genomic analysis indicated that Bra003652 was the candidate gene involved in the swollen root formation in turnips. Our research provides some candidate genes and lays a foundation for the study of the molecular mechanism of swollen root formation.