Four phenotypic Chinese cabbages were used in this experimental study: diploid overlapping type cabbage, diploid outward-curling type cabbage, tetraploid overlapping type cabbage, and tetraploid outward-curling type cabbage—named D2X, S2X, D4X, and S4X, respectively. Diploid overlapping type cabbage and diploid outward-curling type cabbage were produced from diploid high-generation inbred Chinese cabbage (Figures 1A, B). Tetraploid overlapping type cabbage and tetraploid outward-curling type cabbage were produced by crossing common diploid cabbage material BP058, which can produce 2n gametes, with tetraploid cabbage materials Q29 and Q81, which were obtained by mutagenesis (Figures 1C, D). The top of the leaves (the seventh leaf from the inside to the outside) of the four cultivars were taken at the pericardium stage (56 days after planting) (Figure 1E). Three strains were selected, each of which was repeated thrice. All samples were frozen in liquid nitrogen for RNA-seq analysis and qRT-PCR expression analysis.

To a culture dish was added 1 cm² of Chinese cabbage cotyledons. The Chinese cabbage cotyledons were dissociated using 1.5 mL of extraction buffer solution and then cut with a blade. The aforementioned solution was filtered into a 2-mL centrifuge tube and centrifuged for 1000 revolutions for 3 min. After the supernatant was poured out, 50 ul of extraction buffer was added to the precipitate and shaken to suspend the precipitate. Following this, 50–100 ul of dye solution was added to the suspension and dyed for 1 min. A centrifuge tube was put into the flow cytometry (model CytoFLEX) to detect the ploidy of Chinese cabbage, and the diploid variety Chinese cabbage “Jibai 4” was used as the control. Three replicates were taken for each sample.

Total RNA from each leaf sample was isolated using TRIzol reagent (Invitrogen, Carlsbad, CA) according to the manufacturer’s instructions. The cDNA library preparation and transcriptome sequencing followed the related literature (Zhong et al., 2011; Liu et al., 2015). Clean reads were mapped to the cabbage reference genome (brassicadb.org/brad/datasets/pub/Genomes/Brassica_rapa/V3.0/) using TopHat2 software (Kim et al., 2013); only unique mapping reads were retained for calculating gene expression. RNA-seq data analysis was performed according to previously published protocols (Trapnell et al., 2010). DEGs were screened by the following standard: FDR <0.05 and absolute value of log₂ratio ≥1. All genes were aligned against public databases (Nt, Nr, COG, Swiss-Prot, and KEGG) to obtain their putative functions. Similarly, mapping against KEGG genes involved the use of KOBAS v. 2.0 software (Xie et al., 2011) using a blast cutoff e-value of 1 × 10⁻⁵, and the significant enrichment threshold value was padj <0.05 or corrected padj <0.05 for conducting KEGG enrichment analysis. The DEGs were then used to perform co-expression network analysis using the R package WGCNA (Langfelder and Horvath, 2008). All raw sequencing data were deposited at the BioProject under accession code PRJNA944431.

To validate the accuracy of the RNA-seq data, ten differentially expressed genes were randomly selected to perform qRT-PCR, with β-actin as the internal reference gene. The qRT-PCR verifications were conducted as previously described using a Thunderbird SYBR qPCR Mix (Toyobo, Shanghai, China) and LightCycler480II, 384 (Roche). Fragments were 80–150 bp in length. All primers are listed in Supplementary Table S1. Primer Premier 5.0 application software was used to design RT-qPCR primers, which were synthesized by Bioengineering (Shanghai). The gene level was calculated by using the 2^-∆∆Ct method.

All data were presented as means with standard deviation (SD). The data were analyzed using SPSS 17.0 by one-way analysis of variance (ANOVA). Significance statistical analysis was calculated by Duncan’s multiple range test. p < 0.05 indicates a significant difference.

Flow cytometry was used to identify chromosome ploidy for the four cultivars (diploid overlapping type cabbage, diploid outward-curling type cabbage, tetraploid overlapping type cabbage, and tetraploid outward-curling type cabbage). The diploid Chinese cabbage Jibai 4 was used as the control group. The results showed that there was only one main peak in Jibai 4. The fluorescence intensity corresponding to the peak was at the position of 50, and the 2n DNA content was 59.15% (Figure 2A). The fluorescence intensity corresponding to the peak value of diploid overlapping type and outward-curling type cabbage was at the position of 60, and the DNA contents of 2n were 72.18% and 72.72%, respectively (Figures 2B, C). The fluorescence intensity corresponding to the peak value of tetraploid overlapping type and outward-curling type cabbage was at the position of 40, and the DNA contents of 4n were 75.80% and 84.10%, respectively (Figures 2D, E). The diploid Chinese cabbage was identified from tetraploid Chinese cabbage using flow cytometry, which further confirmed the accuracy of the test materials. The primer sequences are listed in Table 1.

To explore the morphological characteristics of four cultivars, the height, width, and weight of the leafy head, the number of leaves, the leaf width and length, and the curvature of the blade tip were systematically measured (Figure 3). The height of the leafy head, number of leaves, and curvature of the blade tip were evidently higher for the outward-curling type cabbage than those for the overlapping type, while the opposite was true for the width of the leafy head, and leaf width and length within the same ploidy.

In the overlapping type cabbage, the height, width, and weight of the leafy head, and the leaf width and length were higher in the tetraploid than in the diploid plants. The height and width of the leafy head, and the number and length of leaves were more pronounced in the tetraploid than in the diploid outward-curling type cabbage.

Transcriptome analysis obtained 114.66 Gb of clean data at an average of 6.92Gb per sample, with Q30 >95.10%. The clean reads of each sample were sequenced with the reference genome, and the alignment efficiency ranged from 82.35% to 91.77%. Gene expression levels for each replicate were assessed using principal component analysis (PCA), which indicated that the four cultivars had clearly differentiated gene expression (Figure 4A). A total of 2,612 and 3,573 DEGs were obtained in S2X vs S4X and D2X vs S2X (Figure 4B), while 28 and 4,151 DEGs were obtained in D2X vs D4X and D4X vs S4X, respectively (Figure 4B). Detailed information regarding RNA-seq data of each sample is listed in Supplementary Table S5.

To further identify the genes closely related to the phenotypes of the four varieties, Venn analysis was performed, which displayed unique and common DEGs in each cultivar (Figure 5A). Furthermore, up- and downregulated genes were also exhibited by Venn diagram. Among these DEGs, only two genes were commonly upregulated, while six were downregulated in the four genotypes (Figures 5B, C).

The Venn diagram showed that 24, 482, and 2,176 genes were exclusively upregulated in the D4X, S2X, and S4X genotypes (Figure 5B), while 25, 619, and 1,369 were exclusively downregulated in the D4X, S2X, and S4X genotypes (Figure 5C), respectively. This suggests that the genotype-specific DEGs might contribute to the phenotypic differences.

Subsequently, all DEGs were subjected to KEGG pathway analysis to identify the major metabolic pathways involved. Upregulated genes were significantly enriched in ribosome, alpha-linolenic acid metabolism, starch and sucrose metabolism, and biotin metabolism (Figure 5D). Downregulated genes were mainly enriched in sulfur metabolism, phytohormone signal transduction, photosynthesis-antenna proteins, stilbenoid, diarylheptanoid and gingerol biosynthesis, limonene and pinene degradation, tryptophan metabolism, and phenylpropanoid biosynthesis (Figure 5E).

Transcription factors (TFs) play a key role in regulating plants’ development and their response to environmental stimuli. The TFs that were identified as DEGs in four different genotypes were highlighted and analyzed. A total of 2,937 TFs were identified, in which the highest rates of TFs belonged to the bHLH, AP2/ERF-ERF, WRKY, MYB, NAC, and C2CH2 families (Figure 6A). Of the 44 identified WRKY DEGs (Figure 6B), BraA04g004020.3C.gene was significantly downregulated only in 4XS, while BraA07g022390.3C. gene, BraA08g017840.3C. gene, and BraA04g028840.3C. gene were significantly upregulated in 2XS and 4XS. For 36 differentially expressed NAC genes (Figure 6C), BraA02g004700.3C.gene, BraA07g030360.3C.gene, and BraA03g003810.3C.gene were highly expressed in 2XD. BraA09g003710.3C. gene, BraA07g034350.3C. gene, BraA05g034120.3C. gene, BraA03g022500.3C. gene, and BraA02g043100.3C. gene were significantly downregulated in 4XD and 4XS. Some 66 MYB DEGs (Figure 6D) were identified, in which more genes such as BraA05g001030.3C. gene, BraA10g000820.3C. gene, BraA08g035740.3C. gene, BraA05g036930.3C. gene, and BraA08g021120.3C. gene were highly expressed in 2XD and 4XD. Meanwhile, 35 C2H2 (Figure 6E), 63 bHLH (Figure 6F), and 48 ERF (Figure 6G) were differentially expressed in four cultivars, indicating that TFs play a significant role in the phenotypic differences among the four cultivars.

In the present study, co-expression networks were built on the basis of pairwise correlations among genes according to the trends of gene expression in all examined samples. The dendogram showed that 34 unique modules were identified, with each module depicted by a different colored branch, and each gene depicted by a leaf (Figure 7A). The gene expression profile of each module was represented by its eigengene—its most notable component. The 34 resulting eigengenes each correlated with unique genotypes due to their genotype-specific expression profiles (Figure 7B).

Notably, five co-expression modules comprised genes that were highly expressed in different cultivars, including MEbrown in tetraploid outward-curling type cabbage (T22, T23, and T24), MEmagenta in tetraploid overlapping type cabbage (T7, T8, and T9), MEpink in diploid overlapping type cabbage (T1, T2, and T3), MEblack in diploid outward-curling type cabbage (T19, T20, and T21), and MEblue in diploid outward-curling type cabbage and tetraploid outward-curling type cabbage. Therefore, each of these five modules identified a specific cultivar or a cluster of genes in two to three similar cultivars. For example, 920 genes involved in the MEbrown module were highly specifically accumulated in tetraploid outward-curling type cabbage, which indicated that this group of genes might be responsible for tetraploid outward-curling head morphotypes. A total of 171 genes involved in MEmagenta were highly specifically accumulated in tetraploid overlapping type cabbage, which indicated that these groups of genes might be responsible for tetraploid overlapping head morphotypes. A total of 195 genes involved in the MEpink module were highly specifically accumulated in diploid overlapping type cabbage, which indicated that this group of genes might be involved in diploid overlapping head morphotypes (Figure 7B). A total of 377 genes involved in MEblack module were highly specifically accumulated in diploid outward-curling type cabbage, which indicated that this group of genes might be involved in diploid outward-curling head morphotypes.

Auxin, acting with other phytohormones, is a key regulator of various developmental processes in plants. Some 113 phytohormone-related DEGs were identified as differentially expressed in four cultivars (Supplementary Table S2). Meanwhile, 62 DEGs common to diploid and tetraploid cabbage (Figure 8A), nine DEGs (three abscisic acids and six auxins) unique to diploid cabbage (Figure 8B), and 27 DEGs unique to tetraploid cabbage were identified (Figure 8C). The DEGs in the auxin hormone pathway, including GH3.17, IAA16, IAA17, LAX1, LAX, and three SAUR21, were evidently highly expressed, while two GH3.10, IAA2, IAA3, LAX3, SAUR20, SAUR23, three SAUR32, and two SAUR50 were lowly expressed in diploid overlapping type cabbage and tetraploid overlapping type cabbage (Figure 8D).

Meanwhile, the abscisic acid-related genes PYL1, PYL4, 2 PYL5, PYR1, SRK2G, AIP1, ABF3, and ABI1 were upregulated, while ABF1, ABF4, PP1CA, PP2CA, PYL1, 2 PYL4, PYL6, 3 PYL8, PYL9, SRK2D, and SR K2E were downregulated in diploid overlapping type cabbage and tetraploid overlapping type cabbage (Figure 8D and Supplementary Table S3). The auxin hormone-related genes IAA16, IAA7, SAUR20, and SAUR24 were upregulated, while IAA12, IAA19, SAUR20, SAUR36, and two SAUR50 were downregulated in overlapping type cabbage compared to outward-curling type cabbage (Figure 8E; Supplementary Table S4). These results indicated that differential expression for phytohormone-related genes may account for phenotypic divergence. Furthermore, to validate the accuracy of the RNA sequencing, qRT-PCR was conducted to examine the levels of ten genes (Figure 9). The results of qRT-PCR accorded well with the expression data obtained by RNA-seq.