Identification of Quantitative Trait Loci for Clubroot Resistance in Brassica oleracea With the Use of Brassica SNP Microarray

Increasing clubroot resistance (CR) of Brassica oleracea by ascertaining the molecular mechanisms has been the key focus in modern B. oleracea breeding. In order to identify the quantitative trait loci (QTLs) associated with CR in B. oleracea, 94 F2 vegetative lines which were developed by tissue culture of selfed seeds from the F1 generation between a clubroot-resistant B. oleracea inbred line and a susceptible line, were identified for disease incidence and six CR-associated traits under a lab inoculation by Plasmodiophora brassicae and were genotyped with the 60K Brassica SNP array. Significant correlations were detected for numbers of fibrous roots and P. brassicae content in roots with disease incidence. Nine linkage groups were constructed from 565 bins which covered around 3,000 SNPs, spanning 1,028 cM of the B. oleracea genome with an average distance of 1.82 cM between adjacent bins. A total of 23 QTLs were identified for disease incidence and the other two correlated traits, individually explaining 6.1–17.8% of the phenotypic variation. Several overlaps were detected among traits, including one three-traits-overlapped locus on linkage group C08 and two important overlapped regions between the two CR-associated traits on C06. The QTLs were compared with known CR loci/genes and the novelty of our QTLs was discussed.


INTRODUCTION
Clubroot, caused by the soil-borne obligate Plasmodiophora brassicae, is a devastating disease in Brassica crops including cabbage (Brassica oleracea L. var. capitata) which is one of the most important vegetable crops in the world (Hirai, 2006;Dixon, 2009). It causes serious yield loss in cabbage since the pathogen always induces galls on the plant root, thus hinders the uptake of water and nutrients and finally leads to abnormal growth (Dixon and Robinson, 2010). It is hard to control clubroot by cultural managements or chemical fungicides due to the long period of the pathogen surviving in soil. Therefore, developing resistance cultivars is the most effective way to control this disease.
In our previous identification, a cabbage inbred line was found to be highly resistant to clubroot disease, thus an F2 segregating population was developed from the cross between this line and a susceptible B. oleracea. On the other hand, a 60K Brassica SNP microarray was successfully released (Clarke et al., 2016) and used for the construction of genetic linkage map with high density of markers in B. oleracea (Mei et al., 2017). In the present study, disease incidence and several CR related traits were investigated in this F2 population in lab tests, and QTLs for these traits were identified by using the linkage group built by data from the SNP microarray. Our study will be helpful for clubroot resistance breeding in B. oleracea.

Plant Materials and Resistance Evaluation
An F2 segregating population, comprising of 94 cloned-lines developed by tissue culture, was developed from hybridization between two inbred lines of B. oleracea, '263' and 'GZ87' with diverse resistance levels against P. brassicae. The development of clones for each F2 genotype was conducted according to Luo et al. (2000). The vegetative plants were transplanted into 90 mm × 80 mm pots after rooting and kept in a climate chamber for 1 week (16/8 h light/dark cycle under 26/20 • C). Then each plant was inoculated by watering 2 mL P. brassicae (the 4 th race) resting spores suspension (2 × 10 8 spores/mL) at the stem base using a pipette method (Carlsson et al., 2004;Nagaoka et al., 2010) and kept in the climate chamber with the same condition as before. Phenotypic data were collected at 6 weeks post-inoculation, including disease incidence (DIC, percentage of diseased plants in total plants), fresh weight per plant (FW), fibrous root weight/root weight (FR/R), length of root (LR), number of fibrous roots (NFRs), ratio of root surface covered with fibrous roots (RFR), and P. brassicae content in roots (PCR).
Four biological replicates of resistance test were conducted in two rounds. Ten to twenty two plants per line were tested in each replicate for the calculation of DIC and five plants for other six indexes. Pearson's simple correlations were calculated between traits of interest via SAS software (SAS Institute, 1999).

Genotyping, Map Construction, and QTL Analysis
Genomic DNA was extracted from young leaf of plant using the CTAB method (Allen et al., 2006). DNA samples of two parental lines and the 94 F2 lines were genotyped by the Brassica 60K Bead Chip Array (Infinium R , Illumina, Inc., San Diego, CA, United States) (Clarke et al., 2016). SNPs were aligned to the reference genomes of B. oleracea (version 1.1) 1 using a local BLAST search. Bins were developed for each chromosome using Perl language by combining SNPs with identical genotypes across the F2 population (Zhang et al., 2014) and then used for genetic map construction in IciMapping version 4.1 with default parameters (Wang et al., 2016). QTLs were analyzed using the inclusive composite interval mapping (ICIM) model in IciMapping with optimized parameters (Step = 1 cM; PIN = 0.005). A permutation test with 1,000 permutations was performed for each trait to calculate the threshold of LOD score at the significance level of P = 0.05.

Phenotypic Performance of Parents and F2 Lines
According to the field identification for three successive years at Fuling, Chongqing, China (E107.6459, N29.5709) where happens serious clubroot disease every year caused mainly by the 4 th race of P. brassicae (identified according to Williams host system, data not shown), 'GZ87' exhibited complete resistance with a disease index of 0, while '263' which is a founder parent in our cabbage breeding program showed stable susceptibility to P. brassicae (with average DIC of 72.6% and disease index of 33.7). Under the lab identification, similarly, 'GZ87' exhibited obviously higher resistance (DIC = 0) than '263' (DIC = 80.8) (Figure 1). Wide variations were detected among the F2 lines for all the seven traits, of which FW, LR, and RFR showed normal distributions, while DIC, FR/R, NFR, and PCR exhibited skew distributions (Table 1). Three traits, particularly PCR and NFR, were found to be significantly correlated with DIC (r = −0.881 and −0.640, respectively) ( Table 2). Therefore, QTLs were subsequently screened for DIC, PCR, and NFR.

SNP Detection and Genetic Map Construction
Among 52,157 SNPs in the Brassica 60K SNP array, 21,646 (41.5%) could be aligned to the B. oleracea genome. These SNPs distributed unevenly on the nine chromosomes of B. oleracea, with the highest number on chromosome C03 (3,558) and the    (Figure 2 and Table 3). The bins showed general coincidence between their genetic and physical positions for all LGs except C01 which has a limited number of bins.

QTLs for CR-Associated Traits
With significant LOD thresholds of 4.4, 4.5, and 4.4 which were set the by 1,000-permutation test in IciMapping, 2, 13, and 8 QTLs for DIC, NFR, and PCR, were detected respectively (Figures 2, 3 and

Comparative Analysis Among CR Loci
According to Bolbase 2 which provides information on syntenic regions between B. oleracea and B. rapa (Yu et al., 2013), a QTL interval on C03 (6.28-6.32 Mb, overlapped by NFR.II-3 and PCR.II-2) in the present study was partially syntenic to 0.59-6.22 Mb on A02 of B. rapa where located the CRc-and Pb(Anju)2linked marker 'm6R' (from 2,112,653 to 2,113,153 bp) (Sakamoto et al., 2008;Nagaoka et al., 2010;Lee et al., 2016). A QTL on C07 (NFR.II-7) was found to be syntenic to 11,077,112-11,414,470 bp on A08 of B. rapa, being partially overlapped with a reported CR loci Crr1 which was fine mapped to 10,692,602-11,617,700 bp on A08 in B. rapa (Hasan and Rahman, 2016). The three important regions on C08 and C06 identified no accordance with reported CR loci.

DISCUSSION
In the past, QTLs were mainly identified based on genetic linkage map constructed by DNA molecular markers such as RAPD, RFLP, AFLP, and SSR markers (Landry et al., 1992;Voorrips et al., 1997;Piao et al., 2002Piao et al., , 2004Suwabe et al., 2006), however, this approach is high-input but low-output in identifying candidate genes. Comparatively, the Brassica SNP microarray which has been widely applied in B. napus could provide high-density genetic maps and thus greatly narrowed the QTL regions for interested traits (Liu et al., 2013;Wei et al., 2016). Since no B. oleracea-specific SNP microarray was available before this study, the 60K Brassica SNP microarray which consisted of probes from both A and C genomes (Clarke et al., 2016) was used in the present 2 http://www.ocri-genomics.org/bolbase/ study. Although nearly 60% (30,511) SNPs were filtered out (could not be aligned to C genome), the 3,218 polymorphic sites from the remaining 40% (21,646) SNPs still enabled us to construct a high-density genetic map of B. oleracea which spanned 1028 cM of the B. oleracea genome with an average of 1.82 cM (0.67 Mb) between neighboring bins. Given the highthroughput and time-saving nature of SNP array, this approach will dramatically accelerate the process of QTL fine-mapping in B. oleracea.
In the studies associated with CR in Brassica, disease index was the most widely used indicator for CR level (Crisp et al., 1989;Landry et al., 1992;Matsumoto et al., 1998;Carlsson et al., 2004;Hirai et al., 2004;Hirai, 2006;Dixon and Robinson, 2010;Hatakeyama et al., 2013;Kato et al., 2013;Liu et al., 2013;Hasan and Rahman, 2016;Lee et al., 2016;Huang et al., 2017;Yu et al., 2017). However, it was hard to judge the disease index of vegetative propagates in the present study since the disease grade of each copy was difficult to class due to the lack of main roots in vegetative plants and the difficulties on distinguishing clubroot tubers from exogenous-hormone-induced calluses in tissue culture. Nevertheless, it was relatively easier to discriminate diseased plants from healthy ones, thus disease incidence was recorded, though its accuracy to indicate the resistance level of host might be not as high as that by using disease index. As supplements of disease incidence, we measured several traits which might be potentially correlated with the resistance of vegetative plants, and we found the P. brassicae content in roots was highly correlated with DIC, being consistent with that in B. rapa (r = 0.95) (Zhu et al., 2015). In addition, the NFR was found to be moderately correlated with DIC and high correlated with PCR. Therefore, the two CR-associated traits particular PCR could be used to indicate the resistance level of plants when needed, for example, in the plants without main roots.
Clubroot-resistant genes or loci were widely identified in B. rapa (Hirai et al., 2004;Piao et al., 2004;Suwabe et al., 2006;Sakamoto et al., 2008;Huang et al., 2017;Yu et al., 2017) and B. oleracea (Landry et al., 1992;Voorrips et al., 1997;Rocherieux et al., 2004;Moriguchi et al., 2010;Nagaoka et al.,   2010; Lee et al., 2016). It was hypothesized that the CR of B. rapa was possibly controlled by qualitative plus quantitative loci (Piao et al., 2002;Suwabe et al., 2006;Sakamoto et al., 2008), while CR trait in B. oleracea might be controlled by QTLs in a quantitative manner (Crisp et al., 1989;Voorrips et al., 1997;Piao et al., 2009;Tomita et al., 2013). In the present study, although only one genetic region on C08 was identified for DIC, many QTLs were identified for the other two CR-associated traits. These QTLs including the two QTLs for DIC explained limited phenotypic variations (6.1-17.8%), suggesting the quantitative nature of these traits and indicated that the pyramiding of these loci may confer durable CR in B. oleracea. Two possible overlaps of our QTLs were detected with previous reported loci after deep comparison of our regions with other studies which provided direct or indirect (primer sequences of linked markers were then blasted to reference genome to get the physical positions) genomic information of CR loci (Matsumoto et al., 1998;Saito et al., 2006;Suwabe et al., 2006;Sakamoto et al., 2008;Nagaoka et al., 2010;Hatakeyama et al., 2013;Kato et al., 2013;Hasan and Rahman, 2016;Lee et al., 2016), however, the QTLs in these two regions only explained 6.7-10.3% of the total phenotypic variation. It is interesting that no overlap was detected for reported CR loci with the three important regions found from C08 and C06 in our study. Although A08 of B. rapa was reported to carry important CR QTLs such as Crr1 and Rcr9 (Hasan and Rahman, 2016;Yu et al., 2017), they are not syntenic to our region on C08. These data suggest that the CR in our resistant B. oleracea GZ87 is possibly controlled by some novel loci. In each of the three genetic regions, we found several pathogen-responsive genes, including genes encoding receptor like proteins (R proteins), receptor binding proteins, auxin efflux transporters, auxin-responsive proteins, RHD3s (root hair defective 3, required for cell expansion and normal root hair development), pathogen-responsive factors (alpha-dioxygenase, calcium-binding, thaumatins, et al.) and transcriptional factors (MYB,ERF,C2H2,et al.). Further study will be carried out to identify the candidate genes and develop linked markers for a pyramiding purpose, since accumulation of QTLs could confer broad-spectrum CR in B. oleracea .