Gene mapping and preliminary functional verification of BcYUCCA6 associated with wrinkled leaf trait in non-heading Chinese cabbage

Non-heading Chinese cabbage (NHCC; Brassica campestris [syn. Brassica rapa] ssp. chinensis) is a cruciferous (Cruciferae) Brassica crop native to China. The phenotypic traits of the leaf, the main agronomic trait-bearing organ of NHCC, are key indicators for assessing the crop's appearance quality^[1]. The molecular regulation of leaf development relies on the synergistic interaction of multiple core regulatory components, which include functional genes, intracellular signaling pathways, and hormone regulatory networks^[2,3]. Previous studies have identified several genes related to leaf wrinkling and curling. In Chinese cabbage, BcpLH participates in the formation of leaf heads by regulating specific subsets of microRNAs, and its downregulation leads to curled and wrinkled leaves^[4]. Overexpression of the sweet potato NAC family transcription factor IbNAC43 induces leaf curling, while knockout significantly alleviates this phenotype^[5].

Plant endogenous hormones are also closely related to the leaf wrinkling phenotype. Relevant studies have shown that mutations and changes in the expression levels of genes involved in plant hormone synthesis pathways or signal transduction pathways could significantly affect leaf morphogenesis^[6]. Auxin is an important plant hormone that acts as a central regulator of leaf growth and development, including leaf initiation, axial growth, and leaf morphogenesis^[7]. Overexpression of Arabidopsis AtYUCCA6 (a key gene for auxin biosynthesis) in potato shows a high-auxin phenotype with narrow and downward-curled leaves^[8]. Mutation of the auxin-responsive protein IAA2 leads to leaf curling in rapeseed^[9]. Gibberellin (GA) also has a significant impact on plant leaf development. Brcl1 is a negative regulator of the GA signaling pathway that inhibits cell proliferation and expansion, and its loss of function causes leaf shape bending in rapeseed^[10].

Bulked Segregant Analysis (BSA) is a classic forward genetic method. Its core idea is to construct extreme trait pools using individuals with extreme phenotypes from genetic segregation populations. Researchers can screen molecular markers closely linked to target traits by comparing allele frequency differences between the two pools. Currently, it has been widely used in trait mapping research of various crops^[11]. In NHCC, this method has been successfully used for gene mapping of traits such as β-caryophyllene synthesis^[12], leaf color change^[13], flower color change^[14], and wax synthesis^[15]. In this study, using the flat-leaf cultivar 'Suzhouqing' (SZQ) and the wrinkled-leaf cultivar 'Xiangqingcai' (XQC) as parents, an F₂ segregation population was constructed. Combined with BSA technology and fine mapping, the key gene BcYUCCA6 for the leaf wrinkling trait was successfully obtained, and its function was initially verified. The results showed that BcYUCCA6 is involved in regulating leaf shape development of NHCC, which lays a foundation for further analyzing the molecular regulatory mechanism of this trait.

In this study, the wrinkled-leaf cultivar 'Xiangqingcai' (XQC) of NHCC was used as the female parent and the flat-leaf cultivar 'Suzhouqing' (SZQ) as the male parent. They were crossed to produce F₁ progeny, which were self-pollinated to generate the F₂ population. The population was planted in fall 2024 at the Baima Research Station of Nanjing Agricultural University in Jiangsu Province for BSA and fine-mapping analysis, with a total of 1,200 plants.

Seeds of 'XQC' and 'SZQ' were germinated for 3 d and then cultured in an artificial climate chamber. At 20 d after germination, leaf samples were taken for paraffin section preparation and hormone content determination. For gene expression pattern analysis, each organ (leaves, heart leaves, hypocotyls, roots, and stems) was collected at 40 d after germination, from three biological replicates.

Phenotype identification and grading standards

After the F₂ population grew for about 85 d, agronomic traits were investigated and statistically analyzed to determine the number of flat-leaf and wrinkled-leaf individuals in the progeny population. According to the surface wrinkle degree, wrinkle depth, and overall undulation characteristics of the leaves of F₂ population plants, the top view phenotype of leaves was divided into five grades, which were used as phenotype identification indicators (Table 1). All experimental data were statistically analyzed using IBM SPSS Statistics 27 software.

BSA analysis for leaf shape trait

After determining the grading of the F₂ segregation population, 30 plants with wrinkle degree grade 1 and 30 plants with grade 5 were selected to form two extreme pools for leaf shape: the flat leaf pool (FL) and the wrinkled leaf pool (WL). Genomic DNA was extracted using the Tsingke DNA extraction kit, and the DNA library was constructed using the TruSeq Library Construction Kit. After library construction, sequencing was performed on the Illumina NovaSeq 6000 platform with a sequencing depth of 30×. The original sequencing data were filtered and then aligned to the reference genome of 'XQC' (Pakchoi-XQC-v1.0, www.tbgr.org.cn) using BWA alignment software (Burrows-Wheeler Aligner, v0.7.13). After alignment, GATK HaplotypeCaller was used to detect SNPs and InDels in each sample^[16]. Finally, SnpEff software was used to annotate all variant sites. Euclidean Distance (ED) was used for candidate interval association analysis, and the square of ED was used for curve fitting in the experiment^[17]. According to the differences in SNP and InDel genotype frequencies between the pools, loci associated with target traits were screened across the whole genome and annotated.

Marker development, genetic map construction, and mapping

In the candidate interval, InDel marker loci that were polymorphic between parents, homozygous in parents, and polymorphic in the population were selected as a set of candidate InDel markers. Sequences of 200 bp upstream and downstream of these loci were selected, and corresponding primers were designed using SnapGene Viewer software (www.snapgene.com) for genotyping of 500 F₂ plants (Supplementary Table S1).

Genomic DNA of leaves from parents F₁ and F₂ plants was extracted using the TPS (1 M Tris-HCl [pH 8.0], 0.5 M EDTA [pH 8.0], 1 M KCl) method. The PCR amplification system included 5 μL Green Taq Mix, 3 μL ddH₂O, 0.5 μL Primer-F, 0.5 μL Primer-R, and 1 μL DNA template. The amplification program was as follows: pre-denaturation at 95 °C for 3 min; 35 cycles of denaturation at 95 °C for 15 s, annealing at the specific Tm value of primers for 15 s, and extension at 72 °C for 10 s; final extension at 72 °C for 5 min. The amplified PCR products were detected by 8% polyacrylamide gel electrophoresis (PAGE) to determine the genotyping results, and the band patterns were counted.

Sequence and expression pattern analysis of candidate genes

The genes in the candidate interval were aligned with the Arabidopsis database (www.arabidopsis.org). TBtools was used to extract the 2,000 bp sequence upstream of the start codon (ATG) of the gene, and the online software PlantCARE (https://bioinformatics.psb.ugent.be/webtools/plantcare/html) was used to analyze the cis-acting elements contained in the promoter region. SnapGene Viewer software was used to align the amino acid sequences. RNA was extracted from the collected samples and reverse-transcribed into cDNA for qRT-PCR analysis. BcActin served as the internal control (Supplementary Table S2). The 2^−ΔΔCᴛ method was used to standardize the relative expression levels of candidate genes^[18]. Each sample was subjected to three biological replicates and three technical replicates. GraphPad Prism 8 was used for data visualization and significance analysis of differences.

Preparation of paraffin sections

Leaves of 'XQC' and 'SZQ' grown for 20 d were fixed in FAA fixative (70% ethanol: 10% formalin: 5% acetic acid: 15% distilled water, v/v/v/v) and stored at 4 °C for 2 d. The fixed samples were paraffin-embedded, and the materials were cut into thin sections (5 μm thick) using a pathological microtome (Leica Instruments Shanghai Co., Ltd., RM2016). The prepared sections were sequentially immersed in environment-friendly dewaxing and clearing solution I for 20 min; environment-friendly dewaxing and clearing solution II for 20 min; absolute ethanol I for 5 min; absolute ethanol II for 5 min; 75% alcohol for 5 min, and rinsed with tap water. The sections were stained using the safranin-fast green double staining method. After staining, the sections were cleared in xylene for 5 min and mounted with neutral gum. After mounting, the leaf tissue structure was observed under a Zeiss upright fluorescence microscope, and ZEISS ZEN3.8 software was used for observation and photography.

Determination of hormone content

A total of 0.5 g of plant leaf tissue powder was weighed, and 5 mL of extract A (ultrapure water : isopropanol : hydrochloric acid = 1:2:0.002, v/v/v) was added. The mixture was treated in a shaker at 100 r·min⁻¹ at 4 °C for 30 min, then 8 mL of dichloromethane was added, and treated in a shaker at 100 r·min⁻¹ at 4 °C for 30 min. The mixture was transferred to centrifuge tubes and centrifuged at 5,000 r·min⁻¹ for 15 min at 4 °C. At this time, the solution was layered, and the powdery tissue was flocculent and located between the two layers of solution. The lower-layer solution was aspirated, dried with nitrogen, dissolved in 1 mL of mobile phase (methanol + phosphate buffer [pH 3.5] = 45:55, v/v), filtered through a 0.45 μm organic phase filter membrane, and stored in a sample bottle at 4 °C for determination using a triple quadrupole liquid chromatography-tandem mass spectrometer (Triple Quad 6500+ Low Mass).

Gene silencing and phenotype observation

According to a previous study by Yu et al.^[19], a 40 bp sequence (5′-TTAACCGACCCCGGTTAGGTCCATTGGAGCTCAAAAATCT-3′) was selected based on the BcYUCCA6 coding sequence, and an 80 bp palindromic sequence was formed by reverse complementation, which was synthesized by Nanjing GenScript Biotechnology Co., Ltd. The synthesized fragment was cloned into the pTY vector to generate pTY-BcYUCCA6. pTY and pTY-BcYUCCA6 plasmids were coated on gold powder, respectively, and 'XQC' with consistent growth for 2 weeks was bombarded with a gene gun to transform the plasmids into the plants. The 'XQC' bombarded by the gene gun was cultured in an artificial climate chamber under the same conditions as described in section "Plant material". The phenotypes of the silenced plants were observed after 2 weeks. RNA was extracted, and the expression level of BcYUCCA6 was detected by qRT-PCR.

The female parent 'XQC' had wrinkled leaves, and the male parent 'SZQ' had flat leaves. The leaf wrinkling trait of F₁ was closer to the female parent 'XQC' (Fig. 1a-c). Longitudinal section observation of leaves by paraffin section showed that the epidermal cells of 'SZQ' leaves were arranged relatively uniformly and neatly. Compared with 'SZQ' leaves, those of 'XQC' exhibited a more locally irregular cell arrangement, particularly in the palisade and spongy tissues (Fig. 1d). Phenotype investigation was carried out on the F₂ segregation population, and 1,200 individuals were divided into five grades according to the grading standards (Fig. 1e). The frequency distribution of leaf wrinkling in the F₂ segregation generation was statistically analyzed, and the normality test of leaf wrinkling trait was performed using IBM SPSS Statistics 27 software. The results showed that the kurtosis was 0.174 and the skewness was 0.194, both < 1, and the coefficient of variation was 33.17% (Supplementary Table S2). The leaf wrinkling phenotype exhibited an approximately normal distribution in the population, indicating that this trait is a complex quantitative characteristic governed by multiple genes (Supplementary Fig. S1).

Figure 1. 

Segregating population construction and related phenotypes. (a) Population construction scheme. (b) Plants (left to right): 'XQC' (♀), F₁ 'SZQ' (♂); scale bar = 10 cm. (c) Leaf comparison of parents and F₁; scale bar = 10 cm. (d) Longitudinal sections of parental leaves (left: 'XQC'; right: 'SZQ'); scale bars = 200 μm (left), 50 μm (right). (e) F₂ population phenotypic grading (Grades 1–5, left to right); scale bar = 10 cm.

BSA analysis and initial mapping of genes controlling the leaf wrinkling trait in NHCC

From 1,200 F₂ population plants, 30 plants with wrinkle degree grade 1 and grade 5, respectively, were selected to form two extreme pools related to the leaf shape trait. The pool sequencing results showed that a total of 195.35 Gb bases were obtained. To ensure the sequencing quality, the original data were filtered to obtain clean reads for subsequent analysis. After filtering, a total of 191.39 Gb clean bases were obtained, with Q30 (Phred score ≥ 30) ranging from 95.15% to 96.28% and GC content ranging from 37.88% to 38.09% (Supplementary Table S3). Both Q20 (Phred score ≥ 20) and Q30 were higher than 95%, and the GC content was stable. The results indicated that the BSA-seq sequencing data were of high quality and could be used for subsequent gene mapping.

Based on these high-quality sequencing data, association analysis was conducted using the Euclidean distance (ED) algorithm based on 5,053,361 high-quality SNP loci identified after filtering. To reduce background noise, the original ED values were squared (ED²), and genomic regions corresponding to SNP loci within the top 1% of ED² values were selected as candidate intervals. A total of 17 regions significantly associated with leaf shape variation were identified, spanning a combined physical interval of 4.08 Mb and containing 753 candidate genes (Fig. 2; Supplementary Table S4). The association map revealed two major regions on chromosome 6 (physical positions: 34.01–35.78 and 46.09–47.17 Mb, spanning 1.77 and 1.08 Mb, respectively), whose ED² values significantly exceeded the predefined threshold, indicating that chromosome 6 harbors key genetic loci governing leaf shape. Although potential association signals were also detected on chromosome 7, their mapping resolution and association strength were comparatively weak. Therefore, chromosome 6 was selected as the focus for subsequent fine mapping.

Figure 2. 

Distribution of ED association values on chromosomes. The abscissa represents the chromosome names. The colored dots represent the ED values of each SNP locus. The black line represents the fitted ED values, and the red dashed line represents the significant association threshold.

Fine mapping of the leaf wrinkling trait gene

Based on the BSA-seq mapping results, we performed fine-mapping verification for both major intervals on chromosome 6 (34.01–35.78 and 46.09–47.17 Mb). Among them, only the 34.01–35.78 Mb interval produced recombinant individuals, while no recombinant plants were detected in the 46.09–47.17 Mb interval with the present population. Therefore, we focused on the interval 34010000–35780000 on chromosome 6 for subsequent fine mapping.

Through 16 pairs of polymorphic InDel markers uniformly developed between parents and F₁, genotyping was performed on 500 F₂ population plants, and 19 recombinant individuals were screened. Combined with the co-segregation analysis of genotypes and phenotypes of recombinant individuals, the target gene was fine-mapped between markers Chr6_34638884 and Chr6_34686067, with a physical distance of approximately 46 kb. This interval contained seven genes (Fig. 3; Table 2). Alignment of candidate gene CDS sequences with the reference genome revealed InDel variations between the parents in Bra06G24347, Bra06G24348, and Bra06G24350. Quantitative expression analysis of the seven candidate genes showed that only Bra06G24348 had a significant difference between the parents, with a significantly higher expression level in 'XQC' than in 'SZQ'. The relative expression levels of other genes were low and had no significant differences (Fig. 4). Based on the above results, Bra06G24348 was identified as a strong candidate gene regulating the leaf wrinkling trait. The homologous gene of this gene in Arabidopsis is AT5G25620 (YUCCA6), so Bra06G24348 was named BcYUCCA6.

Figure 3. 

Results of fine mapping. Lines of different colors represent different genotypes; blue represents 'SZQ', yellow represents F₁, and green represents 'XQC'. The annotation on the left end indicates the parents, F₁, and the serial numbers of F₂ recombinant individuals. The annotation on the right end represents the phenotypes: S is the phenotype of 'SZQ', X is the phenotype of 'XQC'; the black dashed line is the final candidate region.

Figure 4. 

Analysis of gene expression levels within candidate intervals. Values represent the mean ± SD of three biological replicates. *** p < 0.001 (Student's t-test).

Analysis of BcYUCCA6 sequence and expression pattern

To further investigate the differences in BcYUCCA6 between 'XQC' and 'SZQ', the gene was cloned from both genomic DNA and cDNA. The results showed differences at the mRNA level, leading to a change in the C-terminus of the translated protein. Specifically, the 'XQC' protein contains 432 amino acids, whereas the 'SZQ' protein comprises 437 amino acids (Fig.  5a), and this structural difference is not located within the conserved domain (Fig. 5b). Analysis of the promoter region revealed that, compared with 'XQC', 'SZQ' carries a 4 bp deletion (5′-TATA-3′) at position –1,725 bp from the start codon, which precisely removes a TATA cis-element (Fig.  5c). These sequence variations might be related to the leaf wrinkling phenotypic difference between 'XQC' and 'SZQ'.

Figure 5. 

Sequence analysis and organ-specific expression of BcYUCCA6 and leaf endogenous IAA content in two parents of NHCC. (a) Amino acid sequence analysis. (b) Conserved domain analysis. (c) Promoter sequence analysis. (d) Relative expression levels of BcYUCCA6 in different organs. (e) Comparative analysis of endogenous IAA content in leaves. Values represent the mean ± SD of three biological replicates. * p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001 (Student's t-test).

Expression of BcYUCCA6 in 'XQC' and 'SZQ' was determined in different organs using qRT-PCR analysis. The results showed that the expression level of BcYUCCA6 in 'XQC' was higher than that in 'SZQ', with the largest difference in the root organ (12-fold), except in the stem (Fig. 5d). BcYUCCA6 belongs to the YUCCA family of flavin monooxygenase-like proteins, which is a key rate-limiting enzyme catalyzing the conversion of indole-3-pyruvic acid (IPA) to indole-3-acetic acid (IAA) in Arabidopsis^[26]. Therefore, we examined the endogenous IAA content and showed that IAA content in 'XQC' was significantly higher than that in 'SZQ' (Fig. 5e).

Combined with promoter analysis results, it is speculated that the deletion of the TATA element in the BcYUCCA6 promoter of 'SZQ' may reduce transcription efficiency, thereby leading to a significant decrease in gene expression level. In contrast, the intact TATA element in 'XQC' ensures efficient gene transcription, thereby potentially promoting auxin biosynthesis.

Effect of BcYUCCA6 silencing on leaves of NHCC

To further study the function of BcYUCCA6 in regulating leaf shape wrinkling in NHCC, BcYUCCA6-silenced lines were obtained using VIGS. About 15 d after injection of the VIGS vector, the plant leaves faded and showed mosaicism, indicating successful virus infection of the plants. Compared with the control group pTY-Empty, the pTY-BcYUCCA6 plants had weaker blister-like undulation, and the leaves of the silenced plants were flatter than those of the wild-type and control plants (Fig. 6a). RNA was extracted from the leaves of pTY-Empty and gene-silenced plants, and the expression level of BcYUCCA6 was detected by qRT-PCR. In BcYUCCA6-silenced lines #1, #2, and #3, the gene expression level was significantly decreased, with a silencing efficiency of about 87%, indicating successful silencing of BcYUCCA6. Additionally, the expression levels of auxin biosynthesis-related genes IAA1, IAA2, and IAA4 were significantly downregulated (Fig. 6b). To understand the changes at the cellular structure level, safranin-fast green staining was used. By comparing the leaf sections of silenced plants, it was found that the cells of silenced plants were arranged more neatly (Fig. 6c).

Figure 6. 

(a) VIGS silencing of BcYUCCA6 affects leaf development. WT represents wild type, pTY-Empty represents the pTY-Empty plasmid control, and pTY-BcYUCCA6 shows the phenotype of the silenced plants. The scale is 10 cm. (b) Relative expression levels of genes related to auxin biosynthesis in control plants and BcYUCCA6-silenced plants. (c) Leaf sections of wild type, pTY-Empty, and pTY-BcYUCCA6. The scale is 200 μm. Values represent the mean ± SD of three biological replicates. * p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001 (Student's t-test).

Leaf morphology is a pivotal trait of NHCC that directly impacts its agronomic value and market acceptability. Leaf wrinkling, a prominent morphological variation in NHCC, is generally regulated by major genes or polygenic networks, and genetic analysis has verified it as a quantitative trait governed by two pairs of additive-dominant-epistatic major genes^[27]. To elucidate the genetic basis underlying leaf wrinkling, we applied forward genetic analysis. By integrating BSA and fine mapping, we defined the candidate interval as a 46-kb segment on chromosome 6. This interval harbors seven candidate genes. Through qRT-PCR and sequence analysis, we ultimately identified Bra06G24348 as a strong candidate gene, designated as BcYUCCA6.

Notably, BSA-seq detected two major associated intervals on chromosome 6 (34.01–35.78 Mb and 46.09–47.17 Mb) during genome-wide association analysis. We performed fine-mapping verification for both intervals, and only the 34.01–35.78 Mb interval produced recombinant individuals in the present F₂ population, whereas no recombinant plants were detected in the 46.09–47.17 Mb interval with the present F₂ population. In future research, the second major interval and other minor-effect loci will be further studied using enlarged genetic populations (e.g., F_2:3 families) to obtain sufficient recombinant individuals to fully elucidate the polygenic regulatory network of leaf wrinkling in NHCC.

Our previous studies have revealed that genetic mapping of leaf-related traits in NHCC using different recombinant inbred line (RIL) populations identified multiple quantitative trait loci (QTLs) and auxin-related candidate genes. Twenty-seven QTLs controlling 11 leaf-related traits (including petiole length [PL], petiole width [PW], and petiole thickness [PT]) were detected in the 'Suzhouqing' × 'Maertou' population, with auxin-related genes such as SAUR9 and IAA19 annotated in the candidate intervals^[28]; further research on petiole traits using the 'Wutacai' × 'Erqing' RIL population identified 19 stable QTLs and screened BcIAA31 and BcBOP2 as key candidates^[29]. Collectively, these findings have laid a solid foundation for the genetic mapping of leaf traits. In the present study, we adopted BSA-seq to narrow the candidate region associated with leaf wrinkling and subsequently characterized BcYUCCA6 as a strong candidate gene. Our results further enrich the overall understanding of the genetic basis underlying leaf development in NHCC.

Plant hormones are key regulators of leaf shape development, with distinct regulatory pathways mediated by different hormone types^[30]. Brassinosteroids modulate mesophyll tissue arrangement and cell proliferation via transcription factors to influence NHCC leaf wrinkling^[31], while cytokinins govern leaf morphological variation through regulating mesophyll cell division^[32]. Among these, auxin serves as a core hub for leaf development, and our previous QTL studies have confirmed that auxin-related genes (e.g., IAA31, SAUR9) are closely associated with leaf and petiole trait variation^[28,29]. As a member of the flavin monooxygenase family, YUCCA6 functions as a key rate-limiting enzyme in the tryptophan-dependent auxin biosynthesis pathway, catalyzing the conversion of indole-3-pyruvate to indole-3-acetic acid (IAA)^[26]; its expression level is closely correlated with endogenous IAA accumulation. Our results showed that the BcYUCCA6 expression levels and endogenous IAA content were significantly higher in the wrinkled-leaf 'XQC' than in the flat-leaf 'SZQ'. VIGS silencing of BcYUCCA6 further supported its regulatory role for leaf wrinkling (Fig. 6a). The decrease in BcYUCCA6 expression significantly alleviated leaf wrinkling, while the expression of auxin-responsive genes was downregulated (Fig. 6b). The role of the YUC gene family in regulating leaf wrinkling in Brassica rapa was also found by transcriptomic association analysis with other cultivars of NHCC^[33]. In summary, our results suggest that YUCCA6 is a strong candidate gene regulating leaf wrinkling in NHCC.

Furthermore, we find a 4-bp TATA deletion in the BcYUCCA6 promoter of 'SZQ', which is in line with the higher expression of BcYUCCA6 in the wrinkled-leaf of 'XQC' than in the flat-leaf of 'SZQ'. Since the integrity of the TATA box directly affects gene transcription efficiency, this is consistent with previous findings^[34,35]. In addition, differences also exist in the C-terminal sequence of the BcYUCCA6 between the two cultivars. Although this domain is not structurally conserved, further experimental validation is required to elucidate its potential impact on leaf wrinkling.

In summary, this study identified BcYUCCA6 as a strong candidate gene regulating the leaf wrinkling trait, contributing to the research on the molecular mechanism of leaf shape regulation in NHCC, and provided a theoretical basis and gene resources for the molecular breeding application of this trait. Future research will focus on developing molecular markers related to BcYUCCA6, verifying their versatility in different genetic populations, and expecting to achieve rapid distinction between flat-leaf and wrinkled-leaf materials through genotyping.