Phenotypic evaluation showed that the SL varied widely
We measured the SL of 172 RILs from 2016 to 2019, with five replications performed each year, and measured the SL of 520 accessions during 2015 and 608 lines during 2017, with ten replications performed each year. The results showed that SL differs tremendously among rapeseed lines with a 76.26% broad-sense heritability, ranging from 4.81 to 10.89cm in the RIL population and ranging from 3.39 to 12.74cm in the GWAS population, 74.07% of lines concentrated in the range of 5.00–8.00cm (Fig. 1a, Additional file 1: Table S1, Additional file 2: Fig. S1a-b). A correlation analysis showed a strong correlation between SL and thousand seed weight (TSW) and yield per plant (YP), and seed number per silique (SN), seed number per plant (SNPP), siliques per the main inflorescence (SMI), harvest index (HI), and seed oil content (SOC) had a weak correlation with SL (Additional file 3: Fig. S2). The SL variation was further analyzed in Winter, Spring, and Semi-winter subgroups, respectively. The results indicated that there are differences among the three subgroups, but not significant. The SL of Semi-winter cultivars is 6.07 ± 1.23cm, while 5.51 ± 1.45cm and 5.31 ± 1.12cm in Winter and Spring accessions, respectively (Additional file 1: Table S2)
The dynamic observation of silique development was carried out using two extreme materials Z068 (5.45 ± 0.29cm) and Z191 (10.28 ± 1.01cm) who from the RIL population used as representatives of short and long silique, respectively. The results showed that there was no significant difference in SL between long and short siliques in the first 3 days after pollination (DAP). At 9 DAP, there was a significant difference in SL between long and short siliques. The length of the short silique tended to be stable at 18 DAP, but the length of the long silique did not elongate until 27 DAP (Fig. 1b). Further histological observation on the inner and outer epidermis of the long and short siliques of 27 DAP showed that the cell length of the long silique on the outer epidermis was significantly longer than that of the short silique, while the cell width of the long silique was significantly smaller than that of the short silique, and both the cell length and the cell width of the long silique on the inner epidermis were significantly larger than those of the short silique (Fig. 1c-j).
RNA-seq found the plant hormones and the transportation and synthesis of carbohydrates may affect the development of rapeseed silique
After a stringent quality filtering process, 79.37 Gb of clean data were obtained from 12 samples, with a Q30 percentage ≥ 94.80% (Additional file 4: Table S2). The obtained 12 samples of clean reads were mapped to the reference genome sequence of B. napus, and the percentages of mapped reads were similar among the 12 libraries (86.37–90.23%), and 82.71–86.29% of the reads were unique mapped (Additional file 5: Table S3). Based on the mapped results, the FPKM of all genes was counted and found that the log10FPKM value of all samples fluctuated slightly (-2.5–2.5) (Additional file 6: Fig. S3a), and the peak values of most genes were between 0 and 1 (Additional file 6: Fig. S3b). The results of correlation analysis showed that there was a high correlation and good repeatability between biological replicate (Additional file 6: Fig. S3c). After DEGs analysis, 21,482 DEGs were remained for further analysis (the number of DEGs is the sum of each DEG set) (Additional file 6: Fig. S3d-e, Additional file 7: Table S4).
To explore the metabolic pathways enriched for the DEGs, the related DEGs in 12 samples were subjected to KEGG metabolic pathway enrichment analysis. In "T1 vs T2", we found that the DEGs, whether from S-silique or L-silique, were significantly enriched in starch and sucrose metabolism, phenylalanine metabolism, phenylpropanoid biosynthesis, etc. pathways (Additional file 8: Fig. S4). The same enrichment result also appears in "T1 vs T3", and the difference is that the plant hormone signal transduction pathway, which was only enriched in the L-silique in "T1 vs T2", was also enriched in the S-silique in "T1 vs T3". In "T2 vs T3", the DEGs of L-silique were still enriched in plant hormone signal transduction, starch and sucrose metabolism, phenylalanine metabolism, etc (Additional file 8: Fig. S4). In the T3 stage, the origin of DEGs was significantly different from T1 and T2 stages, and the genes related to the cutin, suberine, and wax biosynthesis pathway were significantly up-regulated in this stage. Although there were significant differences in DEGs and enrichment pathways among the three stages, the genes related to starch and sucrose metabolism, plant hormone signal transduction, phenylpropanoid biosynthesis, etc. were differentially expressed in different stages of silique development. These results indicated that the genes related to starch and sucrose metabolism, plant hormone signal transduction, phenylpropanoid biosynthesis, etc. maybe play an important role in silique development.
Co-expression network analysis reveals transcript level differences of photosynthesis and secondary cell wall biosynthesis in long- and short-silique
To identify genes related to SL, we performed a weighted gene co-expression network analysis (WGCNA) using non-redundant DEGs. After using a dynamic tree cutting algorithm, a total of 18 distinct co-expression modules containing 47 to 1787 genes per module were identified, and 1585 uncorrelated genes were assigned into a grey module which was ignored in the following study (Fig. 2a). An analysis of the module-trait relationships revealed that the “lightpink1”(r = 0.97, p = 9e-08), “chocolate3”༈r=-0.86, p = 4e-04༉, “darkgoldnrod4”༈r=-0.75, p = 0.005༉, and “lightblue2”༈r = 0.73, p = 0.006༉ module were highly correlated with the SL in the 12 samples (Fig. 2b,). According to the heatmap of the top 20 genes with the high eigengene connectivity (KME) value in the four modules, the “lightpink1” module had an expression peak in the T2 stage of long silique development, similar to the long silique development pattern, while genes in “chocolate3” module were highly expressed in the T3 stage of short silique development (Additional file 9: Fig. S5a-b). The expression level of top 20 genes in "darkgoldnrod4” module was the highest in the T1 stage of short silique development, while the expression level of “lightblue2” module was the highest in the T3 stage of long silique development (Additional file 9: Fig. S5c-d).
A Gene Ontology (GO) enrichment analysis of the “lightpink1” module genes identified ten significantly enriched GO terms, and most of them related to photosynthesis. Interestingly, photosynthesis-related pathway were also enriched in “lightpink1” module by KEGG pathway enrichment analysis (Fig. 2c, Additional file 10: Table S5). There was no significant enrichment of GO terms in “chocolate3” module. One and 13 GO terms were significantly enriched in "darkgoldnrod4” module and “lightblue2” module, respectively. Among them, most of the enriched terms belonged to “molecular function” and “biological process”, including “monosaccharide transmembrane transporter activity” and “sucrose transport” (Additional file 10: Table S5). Psby (KME = 0.993), encoding a protein in photosystem II, is one of the Hub genes in the “lightpink1” module. In photosystem II, Psby is in close contact with Cytb559, which can protect photosystem II from photoinhibition so that the silique can better for photosynthesis and provide more material and energy for cell proliferation and expansion of the silique pericarp [23]. BnaC09g09210D (KME = 0.982), another Hub gene of the “lightpink1” module, is the homologous gene of AtKNAT7. In Arabidopsis thaliana, AtKNAT7 was a homologousdomain transcription factor of TALE gene family, which is involved in the regulation of secondary cell wall biosynthesis, and its expression is up-regulated by SND1 and MYB46 [24]. The genes with weights value between 0.8 and 1 were screened to construct part co-expression network around hub genes. In the network, multiple genes are involved in cell elongation and expansion, such as DFL2, CSLD3, TCH4, CESA6, etc (Fig. 2d). These results suggest that photosynthesis and secondary cell wall biosynthesis played an important role in the change of silique length during development.
QTL and GWAS Co-located several major loci related to SL
Using linkage mapping, we identified 95 QTLs associated with SL, with a logarithm of the odds (LOD) value above 3.0. Of these, 95 QTLs with phenotypic variation explained (PVE) ranging from 0.02–67.06% were identified on nine chromosomes using four years, BLUE, and BLUP data (Additional file 11: Table S6, Additional file 12: Fig. S6). Three Major QTL were located on A07, A09, and C08, respectively. qSLA7-2 on A07 was a major locus that was stably detected across all environments and explained 56.24–67.06% of the phenotypic variation. qSLA9-1 and qSLA9-3 on A09 also were major loci that were stably detected across all environments but only explained 0.83–5.33% of the phenotypic variation. qSLC8-2 on C08 was another major locus that was stably detected across all environments and explained 9.91–11.31% of the phenotypic variation. In addition, qSLA6-1 and qSLA6-2 on A06 were also detected simultaneously in all environments, but they had a low PVE.
Trait-marker associations were performed using the FarmCPU, Blink, CMLM, GLM, and MLM models in two GWAS mapping populations. Similar results were obtained in the two mapping populations with five models (Fig. 3a-b). In order to facilitate further analysis, combined with Q-Q diagram, we finally choose MLM model as the follow-up analysis model (Additional file 13: Fig. S7a-b). In total, we identified 41 SNPs in 60K population on A06 (1), A07 (6), A09 (15), A10 (1), and C08 (18) chromosomes, whereas 181 SNPs in WRG population were identified on A01 (9), A03 (9), A04 (11), A06 (10), A07 (52), A08 (5), A09 (54), A10 (4), C01 (1), C04 (1), C05 (5), C06 (1), C07 (4), C08 (9), and C09 (6) chromosomes using a threshold of 5% after Bonferroni multiple test correction (Additional file 14: Table S7). The SL was associated with three common significant regions located on A07, A09, and C08 chromosomes, respectively. Among these, 40 SNPs forming a haplotype block on A09 (27.51–28.18Mb) were located in the interval of known QTL for SL [6].
Meta-GWAS and polymorphisms in the candidate region were associated with SL
The meta-analysis detected 85 SNPs associated with SL, of which nine SNPs were undetected in all of the single-population GWAS (Additional file 14: Table S7). 14 and 31 SNPs identified in the spring and semi-winters subgroups were confirmed by meta-analysis, respectively, but only three SNPs identified in the winter subgroup were confirmed by meta-analysis. On the A09 chromosome, the confidence interval of a stable major QTL overlapped with the highly associated region detected by GWAS (Fig. 3c-d). 40 SNPs identified in GWAS were confirmed by meta-analysis, and seven SNPs undetected in GWAS were mined by meta-analysis in this interval(Additional file 14: Table S7). 528 gene symbols were found in this region. Among them, 98 were DEGs in RNA-seq, which listed as candidate genes (Additional file 15: Table S8). LD analysis showed that most peak SNPs were mainly concentrated in block 2 and block 6 (Fig. 3e). The peak SNP (S9_28151819) was involved in a 35-kb LD block (block 6) that encompassed 19 SNP markers, and it was 0.91 kb away from the BnaA9.CP12-2 related to carbohydrate anabolism (Additional file 15: Table S8). In addition, the gene BnaA9.NST2 involved in secondary cell wall biogenesis and Aux/IAA family member BnaA9.IAA30 related to auxin signal were also found in this block. The results of RNA-seq and qRT-PCR showed that the expression of BnaA9.NST2 in short silique was significantly higher than that in long silique, especially in the T2 stage (18 DAP)(Fig. 3f).
In block 2, there are two peak SNPs located inside genes BnaA9.SK21 and BnaA9.TMP-C, respectively. BnaA9.SK21 and BnaA9.TMP-C have the same expression pattern in T1 and T2, but BnaA9.TMP-C has a higher expression level in long silique in T3 (Fig. 3g-h). Two nonsynonymous SNPs, S9_27782829 and S9_27788376, are wthin BnaA9.SK21 and BnaA9.TMP-C, respectively. Further analysis of these two SNPs found that accessions with a AA genotype at the S9_27782829 displayed, on average, 18.56% increase silique length compared with the accessions with a TT genotype. The average silique length of AT genotypes was between the AA and TT genotype. The minor allele (A) was represented in only 17% of the 608 accessions (Additional file 16: Fig. S8a). There was a significant difference between GG and CC genotype at the S9_27788376 (p < 0.01), and the CC genotype will increase silique length (Additional file 16: Fig. S8b).
Similar results were found on chromosome A07, where a stable major QTL Co-located with GWAS. 14 SNPs identified in GWAS were confirmed by meta-analysis, and four SNPs undetected in GWAS were mined by meta-analysis in this interval (Additional file 17: Fig. S9). 102 gene symbols were found in this region, of which 46 were DEGs in RNA-seq, which listed as candidate genes (Additional file 15: Table S8). BnaA7.MYB63, homologous gene MYB63 is a transcriptional regulators specifically activating lignin biosynthetic genes during secondary wall formation in Arabidopsis thaliana [25], was only 0.36 kb away from the significant SNP S7_16015077 with a low expression level during the whole silique development (Fig. 3i). The lignin biosynthesis key gene BnaA7.CCR2 was also detected in this QTL region and BnaA7.CCOAMT, another key gene of lignin biosynthesis, was detected in another highly associated region on A07. Three significant SNPs S7_16214445, S7_16214995, and S7_16215169) are located inside BnaA7.ARF17, and BnaA7.ARF17 is a negative auxin response factor that inhibits downstream auxin-related genes. Further haplotype analysis was focused on the gene BnaA7.ARF17, and four major haplotypes were observed, with low frequency haplotypes (less than five accessions) being omitted (Fig. 4a). We conducted multiple comparison tests of SL, and the results showed that Hap2 and Hap3 had shorter SL than Hap4 (P < 0.05), while Hap1 was an intermediate type (Fig. 4b).
We further analyzed the sequence diversity of the major QTL region on chr A07, A09, and C08 among landraces and pseudo-wild ancestral (European turnip and B. oleracea subspecies) genomes based on previously published data [26] (Fig. 4c-d). The landraces had a lower π value than pseudo-wild ancestral in the major loci region of A09 (27.50–29.40 Mb). On A07, the highest π value of the major QTL region (from approx. 15.80 to 16.4 Mb) was in landraces (1.18 × 10− 3) and pseudo-wild ancestral (1.36 × 10− 3). The sequence diversity of the major loci region on chromosome C08 is consistent with A07 and A09 (Additional file 16: Fig. S8c). These results suggested that these major QTLs might be domesticated and selected during the process of rapeseed domestication from wild type to cultivated rapeseed, resulting in a decrease in sequence diversity.
Lignin biosynthesis plays an important role in silique elongation
The important candidate gene BnaA9.NST2, BnaA7.MYB63, BnaA7.CCR2, and BnaA7.CCOAMT, all of them are lignin biosynthesis related genes, found in the major QTL region of A07 and A09. Meanwhile, the pathway of “Phenylpropanoid biosynthesis” was significantly enriched in different developmental stages of long silique and short silique (Additional file 18: Fig. S10a-c). These results indicate that lignin may play an important role in the development of silique. Thus, we compared the expression levels of genes related to the lignin biosynthesis pathway in different developmental stages of long silique and short silique. An interesting phenomenon is that the expression level of lignin biosynthesis related genes in different stages of the long and short siliques development are significantly different, especially in the T2 stage when the length difference between long and short silique increased sharply, the lignin biosynthesis related genes are highly expressed in short silique (Additional file 18: Fig. S10d).
Similar results were obtained from the observation of the silique pericarp tissue section and the determination of lignin content. The lignified degree increased gradually with the silique developmental process. Especially at the T2 stage, the lignification degree of short silique pericarp was significantly higher than that of long silique pericarp (Fig. 5a-f). Correspondingly, the lignin content of the short silique was significantly higher than that of the long silique at the T1 and T2 stages, especially T2. However, there was no significant difference in lignin content between long and short siliques at the T3 stage (Fig. 5g). Consistent with the results of RNA-seq, the expression levels of the four genes in the long silique were significantly lower than those in the short silique, especially the T2 stage (Fig. 5h). All these results suggest that lignin plays an important role in the formation of long and short silique, especially in the rapid elongation period of the silique.