High-resolution mapping of quantitative trait loci controlling main floral stalk length in Chinese cabbage (Brassica rapa L. ssp. pekinensis)

Background For spring-type Chinese cabbage production, premature bolting refers to the excessive elongation of dwarf stems before harvesting. Although quantitative trait loci (QTL) mapping for bolting-related traits have been studied extensively, the main flower stalk length (MFSL) have been rarely investigated. Two inbred lines, 06–247 and He102, have significant differences in the MFSL. In this study, these two materials were selected as parental lines for the construction of a recombinant inbred line (RIL) mapping population. High-density mapping of QTL for the MFSL was performed based on the deep resequencing of parental lines and specific locus-amplified fragment sequencing (SLAF-Seq) of individual recombination inbred lines. Results An F7 population consisting of 150 lines was developed. Deep resequencing of parental lines produced 21.08 gigabases, whereas SLAF-Seq produced an average of 428.35 million bases for each progeny. The total aligned data from the parental lines identified 1,082,885 high-quality single nucleotide polymorphisms (SNPs) between parental lines. Out of these, 5392 SNP markers with a segregation type of aa×bb and average integrity of > 99% were suitable for the genetic linkage map construction. The final map contained 10 linkage groups (LGs) was 1687.82 cM in length with an average distance of 0.32 cM between adjacent markers. Based on the high-density map, nine QTLs for MFSL were found to be distributed on seven chromosomes, and two major-effect QTLs were identified for the first time. The physical distance between adjacent markers of two major-effect QTLs was 44.37 kbp and 121.91 kbp, respectively. Approximately 2056 and 6769 SNP markers within confidence intervals were identified according to the results of parental line resequencing, which involved 24 and 199 mutant genes. Conclusions The linkage map constructed in this study has the highest density in Chinese cabbage to date. Two major-effect QTLs for MFSL in Chinese cabbage were also identified. Among these, a novel QTL associated with bolting mapped on LG A04 was identified based on MFSL. The results of this study provide an important platform for gene/QTL mapping and marker-assisted selection (MAS) breeding for bolting-resistant Chinese cabbage. Electronic supplementary material The online version of this article (10.1186/s12864-019-5810-2) contains supplementary material, which is available to authorized users.


Background
Heading-type Chinese cabbage (Brassica rapa L. ssp. pekinensis, 2n = 2x = 20), which belongs to the Cruciferae family, is a biennial plant [1]. In China, the area planted with Chinese cabbage accounts for about 15% of the total vegetable cultivated area, ranking first in vegetable cultivated area and yield, respectively [2]. The cultivation system and varieties of Chinese cabbage have been so well developed that year-round cultivation is possible in China. Therefore, the varieties can be divided into three main types according to the culture season, including autumn-type, summer-type and spring-type. However, for spring-type Chinese cabbage cultivation, premature bolting is a persistent problem [3,4], which often makes leaf-head lack commercial characteristics and results in huge economic losses to farmers. In this study, highresolution mapping of quantitative trait loci (QTL) controlling the main floral stalk length (MFSL) in Chinese cabbage would provide a valuable basis for effective strategies to obtain bolting-resistant varieties.
With the development in genomics, marker-assisted selection (MAS) has been used widely in crop breeding for enhancing the efficiency of breeding programs [5]. Using this approach, target traits can be selected indirectly using molecular markers that are closely linked to candidate genes or developed from the gene sequences [6]. However, bolting is a complicated trait and regulated by multiple genes that are largely influenced by environmental cues, such as photoperiod and temperature [7]. Thereby it is difficult to identify the genes or tightly linked marker loci required for MAS breeding. In the advent of molecular marker and QTL analyses, a large number of bolting-related QTLs have been identified in B. rapa [8][9][10][11][12][13][14][15][16][17]. These identified QTLs were based on different segregating populations, in addition, boltingrelated traits were evaluated under different environments and locations. Furthermore, the criterion of bolting-related indices is relatively complex. Among the above-stated bolting-related QTL analyses, most of the indices are related to the flowering character such as bolting time (BT, which is the number of days between sowing and bolting, and bolting is judged by visible flower buds that have emerged in 50% of plants in each line) [11,18], flowering time (FT, which is the number of days from sowing to the first flower opening on each plant) and days-to-flowering (DTF, which is calculated as the day when 50% of the healthy plants had flowered) [14,15,17,19,20]. A few studies used bolting index (BI) combined with different bolting-related factors such as bud emergence and the length of flower stalk as bolting indicators [13,17,21]. Moreover, other investigations used days to bolting (DB) to map the bolting-related QTL, defined as the number of days from seed sowing to 5-cm stalk [16,17]. To date, no uniform international standard for bolting-related indices has been established.
For spring-type Chinese cabbage production, boltingrelated characteristics refer to the elongation of dwarf stem rather than the presence of flower buds. Springtype Chinese cabbage is mature and ready for sale. If the shorter dwarf stem emerged with buds, their commodity value will not be affected, but if longer dwarf stem emerged without buds, their commodity value will not be lost. Thus, the length of dwarf stem is the most important characteristic with regard to bolting-related traits. During reproductive stage, the dwarf stem develops into the main flower stalk. The elongation of dwarf stem to the main flower stalk is the process of bolting, which is greatly influenced by the environment factors such as temperature and photoperiod. In the present study, a B. rapa high-density linkage map were constructed by using specific-locus amplified fragment sequencing (SLAF-Seq) approach of a recombinant inbred line (RIL) population. Since SLAF-Seq is a recently developed low-cost and high-efficient method of highthroughput sequence-based technique, which greatly reduces the complexity of high-quality reference genome libraries [22]. This technology exhibits significant advantages in high-throughput SNP (single nucleotide polymorphism) marker discovery and genotyping for constructing genetic map. In addition, several QTLs for MFSL under two growth conditions were identified. This is the first report for QTLs mapping of MFSL in B. rapa.

Phenotypic analysis of MFSL in a Chinese cabbage RIL population
The phenotypic scores of the two parents and the RIL population were measured (Additional file 1: Table S1). The phenotypic performance of the RIL population was continuously distributed and relatively consistent during the two years, suggesting that the trait was inherited in a quantitative manner. The phenotypic data of the RIL population was highly variable, ranging from 3.00 cm to 55.50 cm and 4.10 cm to 66.30 cm in 2016 and 2017, respectively. While MFSL between year 2016 and year 2017 was highly correlated (r = 0.90). The results of frequency distribution analysis of MFSL are shown in Fig. 1. A comparison of the parental lines showed significant differences (p < 0.01) within two years of cultivation ( Fig. 1).

High-throughput sequencing and genotyping of the RIL population
For the two parental lines, He102 and 06-247, the average sequencing read number was 70.40 million, representing 21.08 gigabases (Gb). The Q30 ratio was 90.65%, and average GC content was 38.30%. The average mapped reads and genome coverage for the two parental lines was 97.12 and 86.63%, respectively. However, for each progeny, the average read number was 2.14 million, representing 428.35 million bases. The Q30 ratio was 95.06%, and average GC content was 40.98% for the RI line ( Table 1). The total number of aligned reads from the parental lines identified a total of 1,082,885 highquality SNPs between parental lines, which included 805,813 SNPs that matched the aa×bb segregation pattern. After filtering out lower-quality SNP markers, a total of 5392 high-quality SNP markers were used for linkage map construction. The average sequencing depth for the markers was 33.73× in the parental lines and 27.98× in the offspring, and the average integrity was 99.96%.

Construction of a high-density linkage map
A total of 5392 markers were assigned to 10 linkage groups (LGs) (Fig. 2). The map spanned a total of 1687.82 cM, with an average distance of 0.32 cM or 49 kbp between adjacent markers ( Table 2). The largest LG was LG A09 with 811 markers and a total length of 225.79 cM and an average distance of 0.28 cM between adjacent markers. The smallest LG was LG A07 with 294 markers and a length of 89.76 cM. On the other hand, the marker density could be deduced based on physical distance. For that matter, the densest LG was LG A08, with average physical distance of 36.62 kbp between adjacent markers, while the least dense LG was LG A07, with average physical distance of 88.01 kbp between adjacent markers. The degree of linkage between markers was reflected by a gap ≤5 cM ranging from 99.81 to 100%, with an average value of 99.96%. The largest gap on the map was 7.06 cM, located within LG A02.
Analysis of the segregation of all 5392 mapped loci showed that 5239 (98.16%) significantly deviated (P ≤ 0.05) from the expected 1:1 Mendelian segregation ratio. Only 99 markers showed significant (0.05 < P < 0.001) segregation distortion. These segregation distortion markers (SDMs) were located unevenly within the 10 LGs (Table 2). LG A05 had the highest number of SDMs at 31, and LG A08 had no SDMs.
Correlation analysis between the genetic location and the corresponding physical position of mapped SNP markers is an important indicator for evaluating genetic maps. Collinearity analysis revealed that > 90% of the SNP loci on the linkage map were in the same order as those on the corresponding chromosomes of the  physical map (Fig. 3), making the annotation of genes within QTL intervals feasible.

Mapping of MFSL
Based on the high-density genetic map, nine QTLs associated with MFSL were mapped on seven LGs through inclusive composite interval mapping (ICIM), and QTLs were not detected on LG A05, 07 and 08 (  For spring-type Chinese cabbage production, premature bolting has remained a major challenge. It often causes excessive elongation of the dwarf stem and results in the loss of commercial value for spring-type Chinese cabbage leaf-head. Because the elongated dwarf stem would develop into main flower stalk, we selected parental lines with significant differences in MFSL (Fig. 1) for mapping population construction. However, previous studies about QTL mapping of bolting-resistant traits seldom focused on MFSL.
To construct high-resolution genetic map, a total of 150 RILs were used in our study. To obtain high-quality SNP markers for linkage map construction, the method of combining deep resequencing on parental lines and SLAF-Seq on offspring individual line was conducted.
Here, it was unnecessary to sequence the parental lines for linkage map construction and QTL mapping. However, more polymorphic markers were required for fine mapping of major-effect QTL and map-based cloning of candidate genes. In fact, our strategy resulted in 1,082, 885 SNP markers between parental lines, of which 805, 813 showed an aa×bb segregation pattern. These highquality SNP markers provide a marker pool for fine mapping and map-based cloning of candidate genes in subsequent studies.
To construct a reliable linkage map, the isolated aa×bb pattern SNP markers still need stringent screening. High-depth detection in parental lines could ensure accurate genotyping of offspring individuals. In this research, the parameter setting on the depth of the selected SNP markers was not only larger than 14× in parental lines, but also larger than 5× in at least 90% (the value of integrity) of offspring lines. Actually, the average depth of the selected markers in two parental lines and offspring was 34.18× and 27.98×, respectively. It was significantly larger than previous study provided by Yu et al. [2] in which the average depth of parental lines and offsprings was only 21.96× and 10×, respectively. The integrity of all selected markers was 99.96% on average, significantly > 80% in the previous one [2]. Otherwise, to make sure that the annotation of genes within QTL intervals was feasible, collinearity analysis was performed (Fig. 3). Finally, the linkage map included 5392 high-fidelity SNP markers on 10 LGs, with total map length of 1687.82 cM and an average distance of 0.32 cM between adjacent markers ( Fig. 2 and Table 2), which were significantly larger and denser than previously constructed bin map [2]. Significant differences in total map length were existed between the present study and Yu et al. [2]. It might be caused by the difference in parental lines, mapping population, population size and total marker number [23,24]. The parental lines we selected showed great polymorphism, the RIL population could produce a large number of recombination events, the size RIL150 mapping population was larger than previous DH100, and the total number of mapping marker in present study was also significantly more than previously one [2]. To our knowledge, this map has the highest marker-density to date among individual experimental Chinese cabbage genetic maps.
Based on the present high-density linkage map, highresolution mapping on QTLs of MFSL was performed. In a previous study, QTL mapping of bolting-related traits in B. rapa showed that all 10 LGs cover QTLs, except for LG A04 [17]. This suggested that a novel QTL related to bolting was identified based on present phenotypic indices, namely, the length of main flower stalk. MFSL10 was the other major-effect QTL found in the current study, it might be localized to the bottom of LG A10 according to the ratio of peak to total map length (115/128), which may be applicable to the qDB10-2 previously determined [17]. The flower stalks here were similar to the inflorescence stems in Arabidopsis. In previous studies involving QTL mapping of boltingresistant traits, the indices such as BT [11,18], FT [14,15,17,19,20], DB [8,17], and BI [13,17,21] were concerning whether with flower buds or first flower opening. However, extensive investigations on Arabidopsis revealed that the development of flower and inflorescence stem involves different molecular mechanism [25], and both processes are influenced by the same environmental factors such as temperature and illumination. QTL mapping for plant height (PH) in B. rapa was conducted using chromosome segment substitution lines (CSSLs), and a total of 14 QTLs for PH were identified, but none were located on chromosome A04 [16]. In the study, the PH was measured from the ground to the top of the stem on the day the first flower blossomed. The measuring method of MFSL was similar to PH, but the measuring time was different. Because the first bud blossomed at different times for different lines, the PH value for different lines could not be compared with each other. The notion of PH was related to the blossom of first flower, it might be closely linked with flowering time. The annual harvest time of spring-type Chinese cabbage is fixed. One of most important commercial character for harvested spring-type Chinese cabbage was whether the dwarf stem was overly elongated or not. The excessive elongation of dwarf stem would cause a decrease in the product value. In practice, mature spring-type Chinese cabbage leaf-heads with shorter dwarf stems were of higher demand than those with longer dwarf stems. Therefore, it is necessary to identify markers closely linked with shorter dwarf stems, which will develop into the main flower stalk, for breeding bolting-resistant spring-type Chinese cabbage by MAS. The present work focused on achieving this goal. Moreover the confidence intervals of two stable QTLs identified in present study were away from distance of SDM, it made the results more credible.
The regulatory mechanism of flower stalk elongation is relatively complex. To understand the molecular mechanism underlying flower stalk elongation, candidate genes should be isolated in the future. Candidate genes were screened between the confidence intervals of major-effect QTLs, namely MFSL4 and MFSL10. SNP markers and mutant genes between the confidence intervals were identified. In our future study, we plan to conduct fine-mapping and map-based candidate gene identification.

Conclusions
In this study, a high-density genetic map of B. rapa was constructed using the SLAF-Seq technique. To our knowledge, this map has the highest marker density in B. rapa to date. Furthermore, nine QTLs associated with MFSL, including two major-effect QTLs, were identified. Notably, MFSL04 is a novel QTL associated with bolting mapped on LG A04 based on the MFSL. SNP markers and mutant genes between the confidence intervals were also analyzed. These results lay a foundation for finemapping and map-based candidate genes identification for flower stalk length.

Mapping population and phenotyping
A RIL population consisted of one hundred and fifty F 2:7 recombination inbred lines derived by single-seed descent from a cross between two parents, 06-247 and He102, was used. 06-247, the maternal parent with shorter dwarf stems, was selected from the selfing progenies of Japanese hybrid JianChun (JingShenCai No. 2004007, imported from TAKII Seed Company, Japan). He102, the paternal line with longer dwarf stems, is a selfing progeny of a local variety, Henan Erbaotou, collected from Henan Province in China [26].
Phenotypic experiments were conducted at the Core Experiment Station of the Vegetable and Flower Institute of Shandong Academy of Agricultural Sciences, Jinan, China (36.67 N, 116.98 E) during two spring growing seasons from 2016 to 2017. The seeds of each line were germinated at 28°C for 24-48 h in darkness, and the germinated seeds were sown into 7-cm pots on January 3, 2016. Five plants were selected as one independent biological replicate, and three biological replicates were used for randomized complete block design. The pots were directly placed in a solar greenhouse. In 2016, the greenhouse was covered by a layer of thermal insulating layer when the natural temperature descended below 0°C , at night, whereas in 2017, the greenhouse was covered by two layers of plastic film that were separated by a 50-cm interval to keep warm. Generally, the condition in the greenhouse was 6-8 h day (highest to 20°C) /16-18 h night (lowest to 3°C) with a light intensity of 8000-16,000 lx during both growing seasons before transplanting. On March 10th of each year, all plantlets were transplanted in a pest-prevent net (φ0.8) room until seeds were harvested. The main flower stalk was measured from the ground to the top of the stem after 30 days of transplanting at each growing season. The mean phenotypic data, standard deviation and correlation analysis were performed by microsoft excel. The mean phenotypic data for three replicates of the individual was used for further analysis. The lines with phenotypic data less than three replications were deleted and a final 136 lines were used for QTL mapping.

DNA extraction, library preparation, and sequencing
For DNA extraction, the leaves of the parents and the RIL population were sampled and frozen in liquid nitrogen. High-quality genomic DNA was extracted using the DNeasy Plant MaxiKit (Qiagen, USA) and stored at − 20°C prior to preparation of a sequencing library.
Two different sequencing library construction strategies were used for parent lines and offsprings. For parent lines, standard paired-end (PE) libraries with 350-bp inserts were constructed following the manufacturer's instructions (Illumina, USA). Briefly, extracted gDNA was sheared ultrasonically using a Covaris system (Applied Biosystems, USA). Then, the fragmented DNA was purified with the QIAquick PCR Purification Kit (Qiagen, USA), followed by end-repair and index adapter ligation using the TruSeq PE Cluster Kit v3-cBot-HS (Illumina, USA). The resulting libraries were sizeselected on a 2% agarose gel. The fragments between 300 bp and 500 bp were collected with the QIAquick Gel Extraction Kit (Qiagen, Germany) and used for amplification by polymerase chain reaction (PCR). The resulting PCR products were purified further and used for the sequencing library. For the RIL population, the SLAF-Seq strategy with minor modifications was used according to Sun et al. [22], in which two enzymes (RsaI and HaeIII, New England Biolabs, NEB, USA) were used to digest the genomic DNA. The details of the following operation were referred to in Yu et al. [2]. All highthroughput sequencing was performed on an Illumina HiSeq 2500 system (Illumina, USA) according to the manufacturer's recommendations. The sequencing and the following bioinformatics analysis were conducted by Biomarker Technologies (Beijing, China).

Sequence data analysis
The sequence data of the parental lines and each RIL were distinguished by Illumina Casava 1.8 in FASTQ format. The raw reads were first sorted according to barcodes. Then a series of quality control (QC) procedures were adopted. QC standards were performed using the Trimmomatic program as follows [27]: (1) Reads with > 10 nt aligned to the adapter were removed, (2) Reads containing the RsaI or HaeIII sequences were removed, (3) Reads with unidentified nucleotides (Ns) > 10% were removed, (4) Reads with > 50% of bases having a Q phred < 20 were removed (Q phred can be used to represent average error rate (e), the formula between Q phred and e is as follows: Q phred = − 10*log 10 e.). The sequencing depth for each base and total reads coverage compared to the reference genome were calculated based on the alignments.

Mapping reads and SNP calling
Sequence reads were aligned to the B. rapa reference genome V1.5 (http://brassicadb.org/brad/) [28] with the Burrows-Wheeler Aligner (BWA) using default parameters. Briefly, unmapped reads or the reads mapped to multiple locations were excluded. The aligned reads considered to be PCR duplicates were removed using the MarkDuplicates in the Picard software package. At last, a series of processes were performed using the Genome Analysis Toolkit (GATK) to identify polymorphic sites: (1) Regions near short insertion or deletion (InDel) were realigned locally with IndelRealigner; (2) Base quality scores were recalibrated with Base Recalibration; (3) Variants such as SNPs and InDels were extracted with UnifiedGenotyper; (4) SNP clusters (two SNP in between a 5-bp region) and SNPs within 5-bp regions adjacent to In/Dels were removed. The resulting SNP markers were used for the following genotyping and mapping analysis.

SNP marker selection and genotyping of RIL population
To obtain high-quality SNP markers for genetic map construction, all SNPs identified among parental lines and the offspring were filtered as follows: markers with three or more kinds of base type, markers showing no polymorphisms between parental lines, markers not detected or with a depth of < 4× in the parental lines. The remaining polymorphic markers were potential markers suitable for subsequent analysis and they were classified into eight segregation patterns (ab×cd, ef × eg, hk × hk, lm × ll, nn × np, aa×bb, ab×cc, and cc × ab). The mapping population was a RIL population, so tags showing the aa×bb segregation pattern were used for genetic map construction. For high-confidence genetic map construction, the resulting markers with integrity (defined as the percentage of determined markers for all mapped markers in a given line) < 90% and sequencing depth less than 14× in the parental lines as well as < 5× in the offspring were further filtered.

High-density genetic map construction and segregation distortion markers
Based on the genotyping data of 150 RILs, a highdensity genetic map comprising 10 LGs was constructed by a newly developed HighMap strategy [29]. The detailed protocol of HighMap contains four modules, designed for linkage grouping, marker ordering, error genotyping correction and map evaluation, respectively. In order to produce a high-confidence genetic map, the ordering module and error genotyping correction module were conducted iteratively. Map distances were estimated using the Kosambi mapping function [30]. Segregation distortion markers were unavoidable. Marker segregation ratios were calculated using the chisquare test. The markers showing significant (p < 0.05) segregation distortion were initially eliminated when the linkage map was constructed. Error genotypes of segregation distortion markers were corrected further by the SMOOTH software [31] according to parental contribution of genotypes, and then missing genotypes were imputed based on a k-nearest neighbor algorithm [32]. They were then added as accessory markers.

Collinearity analysis between the genetic and physical maps
All the sequences of SNP markers that were constructed in the linkage map were aligned back to the physical sequence of the B. rapa genome through local Basic Local Alignment Search Tool (BLAST) to confirm their physical positions in the genome. Software CIRCOS 0.66 was used to compare the collinearity of markers based on their genetic positions and physical positions [29].

MFSL evaluation and QTL mapping analysis
The MFSLs of the two parents were compared by the Student's t-test at the 5 and 1% levels of probability. The frequency distribution of MFSL in RIL population was analyzed. Inclusive composite interval mapping in the bi-parental populations (BIP) model of QTL IciMapping software v4.0 [33] was used to detect the additive QTLs for MFSL with the P values for entering variables (PIN) = 0.001. Because many adjacent markers in the SNP linkage map had a map distance of < 1 cM, the scanning step was 1.0 cM. LOD significance thresholds for QTL peaks were determined using 1000 permutations. The percentage of variance explained and the additive effect for each QTL were also estimated.

Additional files
Additional file 1: Table S1. Phenotypic analysis of MFSL in two parents and RIL population of Chinese cabbage. (XLSX 95 kb) Additional file 2: Table S2