In our previous study (Gao et al., 2020), 14 non-heading mutant plants were harvested in an EMS-induced mutagenic population. For this study, we selected two mutants with extremely similar phenotypes, nhm4-1 and nhm4-2.

For genetic analysis, the mutants were crossed with the wild-type “FT” to obtain F₁, F₂, and BC₁ populations. To investigate the genetic characteristics of the mutants, we recorded the phenotypes of each plant in each generation and analyzed the separation rate of the F₂ and BC₁ populations using the Chi-square (χ²) test.

To detect the allelism of the two mutant genes, we conducted an allelism test. Mutants nhm4-1 and nhm4-2 were used as parents for hybridization, and the phenotypes of their hybrid progeny were observed and recorded.

A modified MutMap method (Abe et al., 2012) was applied for fine mapping and identification of candidate genes for nhm4-1. For MutMap, the DNA of 50 mutant plants in the F₂ population was mixed equally as the mutant pool. DNA from the two parental plants and the mutant pool asextracted from fresh leaves at the rosette stage using a DNA secure plant kit (Tiangen, Beijing, China) and resequenced with a NovaSeq 6,000 sequencer (Illumina, San Diego, CA, United States of America).

Low-quality data were filtered from the raw data according to the filtering criteria of our previous study (Gao et al., 2020) to obtain clean reads. Clean reads were mapped to reference genome sequences using BWA software (Li and Durbin, 2010), and SAMtools (Li et al., 2009) was used to sort the alignment file. Insertions and deletions (INDELs) and single-nucleotide polymorphisms (SNPs) were identified using GATK software (McKenna et al., 2010) and ANNOVAR software (Wang et al., 2010). Circos software (Krzywinski et al., 2009) was used to map variation information in the genome. The ΔSNP index across the chromosomes of the B. rapa genome was obtained using sliding-window analysis (with a five-SNP window size and one SNP for each step).

Kompetitive Allele-Specific PCR was developed for the genotypic assay to detect the co-segregation of each SNP and to confirm the nhm4-1 candidate gene. The allele-specific primers used are shown in Supplementary Table S1. A total of 184 F₂ plants were used for KASP genotyping. Of these, 48 plants showed a mutant phenotype and 136 exhibited the wild-type phenotype. The assay was carried out at the Vegetable Research Center of the Agriculture and Forestry Academy in Beijing.

The coding sequences of candidate genes were amplified with specific primers (Supplementary Table S2) in the wild-type “FT,” nhm4-1 mutant, and nhm4-2 mutant plants. Primer 5.0 was used to design specific primers and perform gene cloning following the methods of Gao et al. (2020). Sequencing was performed using the Sanger method at GENEWIZ (Suzhou, China). The sequences were aligned using DNAMAN V6 software (Lynnon BioSoft, Montreal, QC, Canada).

Ent-Copalyl diphosphate synthase 1 activity in the leaves from the wild-type and nhm4-1 mutant plants was evaluated using the Plant CPS enzyme-linked immunosorbent assay (ELISA) Kit (Meimian Biotech Co., Ltd., Jiangsu, China) via a double antibody sandwich method following the manufacturer’s instructions. Each material was performed for three biological repeats, and three times the technical repeats were performed in each biological repeat.

For analysis of the relative expression levels of the candidate gene, total RNA was extracted from the cotyledons, first true leaves, third true leaves, sixth true leaves, rosette leaves, and heading leaves of wild-type “FT,” mutant nhm4-1, and mutant nhm4-2 plants using an RNA extraction kit (Aidlab, Beijing, China). First-strand cDNA was synthesized using the FastQuant RT kit (Tiangen, Beijing, China) and quantitative real-time PCR (qRT-PCR) was performed using an Ultras SYBR Mixture (CWBIO, Beijing, China) and the Quant Studio 6 Flex Real-Time PCR System (ABI, Los Angeles, CA, United States of America). The Actin gene (F: 5′-ATCTACGAGGGTTATGCT-3′; R: 5′-CCACTGAGGACGATGTTT-3′) was used as the reference gene. Each experiment was independently performed with three technical replicates and three biological replicates. The relative gene expression levels were calculated using the 2^−ΔΔCt method (Livak and Schmittgen, 2001). The primer sequences used for qRT-PCR amplification are shown in Supplementary Table S3.

According to our previous methods (Gao et al., 2020), the content of endogenous GA in the leaves of wild-type “FT,” mutant nhm4-1, and mutant nhm4-2 plants was determined by liquid chromatography-tandem mass chromatography (Chen et al., 2012). The response of the mutants nhm4-1 and nhm4-2 to GA was determined by spraying exogenous GA₃ solution (500 mg/L). The plants were sprayed with GA₃ at 3-day intervals once the cotyledons were fully expanded, and the treatments ended before the rosette stage. Seedlings sprayed with an equal volume of ddH₂O without GA₃ were used as controls.

The phenotype of the mutant nhm4-1 was highly consistent with that of the mutant nhm4-2. The leaves of the mutant plants showed geotropic growth throughout the developmental stages and could not form leafy heads at the heading stage, unlike those of the wild-type “FT” plants (Figure 1).

As shown in Table 1, the phenotype of all plants in the F₁ generation was consistent with that of the wild-type “FT.” The segregation ratio of the F₂ generation was 3:1, while that of the BC₁ generation was approximately 1:1. These results indicate that the mutation phenotype of the nhm4-1 and nhm4-2 mutants was controlled by a single recessive nuclear gene.

We crossed mutants nhm4-1 and nhm4-2 to detect the allelism of the two. Both F₁ populations exhibited a mutant phenotype after hybridization, indicating that the mutant genes of nhm4-1 and nhm4-2 are allelic (Figure 2).

The wild-type “FT,” the mutant nhm4-1, and a mutant pool containing 50 homozygous recessive F₂ mutant plants were resequenced, resulting in 94,719,968, 118,123,996, and 270,862,892 clean reads, respectively. A total of 98.00, 98.43, and 98.91% of the clean reads in the wild-type “FT,” mutant nhm4-1, and mutant pools, respectively, were mapped to the Chinese cabbage v. 3.0 reference genome (BRAD¹). Based on the alignment to the reference genome sequence, the mutation analysis software GATK (McKenna et al., 2010) was used to extract all potential polymorphic SNP sites in the genome. Circos software was used to draw the variation information on the genome, and it was found that SNP was mainly distributed on chromosome A03 and A09 (Supplementary Figure S1). This was followed by further filtering and screening, after which 1,587 high-quality SNPs were obtained.

When the SNP index was 0.95 as the threshold, we located a 1.59 Mb (867,020-2,457,084) candidate region on chromosome A09 (Figure 3). Five SNP mutations occurred in the exon, of which only two SNPs (SNP A09, 900,112 and SNP A09, 1,723,490) caused non-synonymous amino acid changes (Table 2). SNP A09, 900,112 (C–T) was located in the BraA09g001440.3C and SNP A09, 1,723,490 (G–A) was located in the BraA09g002790.3C.

To confirm the candidate SNP, primers were designed based on the mutation information of these two SNPs and then applied to F₂ populations. Genotyping analysis was performed using KASP, and the association between these two SNPs and the mutant phenotypes was verified. The genotypic assay showed that SNP A09, 900,112 of BraA09g001440.3C was the T:T genotype in the 48 mutant phenotypic plants and the C:T or C:C genotypes in the 136 wild-type phenotypic plants (Supplementary Table S4), implying that SNP A09, 900,112 co-segregated with the mutant phenotype. However, a recombinant was found at SNP A09, 1,723,490 of BraA09g002790.3C. The A:A genotype and the A:G genotype were detected in the mutant phenotypic plants. Thus, this SNP did not co-segregate with the mutant phenotype. These results confirm that BraA09g001440.3C harbors SNP A09, 900,112, and is the most likely candidate gene of the nhm4-1 mutant. Gene annotation indicates that BraA09g001440.3C is a homologous gene of Arabidopsis CPS1 (At4g02780) and encodes the CPS1 enzyme, which catalyzes the conversion of geranylgeranyl pyrophosphate (GGPP) to copalyl pyrophosphate (CPP) in GA biosynthesis. In this study, the candidate gene of the nhm4-1 mutant is referred to as BrCPS1.

Since MutMap and KASP analyses supported BrCPS1 as the most likely candidate gene of Brnhm4-1, we cloned the coding sequence of BrCPS1 from the wild-type “FT,” mutant nhm4-1, and mutant nhm4-2. Gene annotation showed that BrCPS1 was 7,928 bp in length and consisted of 15 exons (Figure 4). Sequence alignment of the cDNAs and deduced amino acids between the wild-type and mutants is shown in Figure 5. Sequence comparison showed that a single base substitution occurred at position A09, 900,112 (C to T) in nhm4-1, resulting in an amino acid to change from leucine (L) to phenylalanine (F). However, the BrCPS1 clone in mutant nhm4-2 plants showed a single-nucleotide mutation in the fourth exon (A09, 898,783; G–A), resulting in an amino acid to change from glycine (G) to aspartic acid (D). This differs from the mutation site of the nhm4-1 mutant, suggesting that nhm4-1 and nhm4-2 are allelic mutations of the same BrCPS1 gene.

In this study, we also cloned the candidate gene BraA09g002790.3C (SNP A09, 1,723,490) in mutants nhm4-1 and nhm4-2. There was a single base mutation (G–A) in nhm4-1 at position A09, 1,723,490, which is consistent with the MutMap results. There was no change in the mutant nhm4-2. These results further confirm that BrCPS1 is the gene responsible for the non-heading phenotype.

To investigate whether the activity of CPS1 was changed in mutant nhm4-1, we measured the activity of CPS1 in wild-type and mutant nhm4-1. The activity of CPS1 was significantly decreased in the mutant nhm4-1 compared with the wild-type (Supplementary Figure S2).

To investigate whether the mutation site affects gene expression, we used qRT-PCR to analyze the expression level of BrCPS1 in the cotyledons, true leaves, rosette leaves, and heading leaves of wild-type “FT,” mutant nhm4-1, and mutant nhm4-2 plants. The expression levels of BrCPS1 were reduced at all stages of leaf development in both mutants, especially in the rosette leaves, which differ from wild-type “FT” results (Figure 6).

We determined the GA content of 18 endogenous GAs (GA₁, GA₃, GA₄, GA₅, GA₆, GA₇, GA₈, GA₉, GA₁₂, GA₁₅, GA₁₉, GA₂₀, GA₂₄, GA₂₉, GA₃₄, GA₄₄, GA₅₁, and GA₅₃) in the leaves of wild-type and mutant plants. GA₅, GA₆, GA₇, GA₄₄, and GA₅₃ were not detected in either the wild-type or mutant strains (Figure 7). GA₁₂ was detected only in the wild-type and not in either of the mutants. The levels of GAs (GA₉, GA₁₅, GA₂₀, and GA₂₄) in the bioactive GA biosynthesis pathway were significantly decreased in the mutants than in the wild-type.

Based on the determination results of endogenous GA content, we investigated the responses of the mutants to exogenous GA₃ application. After exogenous spraying of GA₃, we observed that the leaves of the mutants had grown upward, similar to the wild-type leaves in the rosette stage, which demonstrates a tendency to form a leafy head (Figure 8). Consequently, leafy head formation is related to GA biosynthesis or the deactivation pathway.