Key genes and mechanisms underlying natural variation of silique length in oilseed rape (Brassica napus L.) germplasm

Quaid Hussain,Jiepeng Zhan,Huabing Liang,Xinfa Wang,Guihua Liu,Jiaqin Shi*,Hanzhong Wang*

Oil Crops Research Institute of the Chinese Academy of Agricultural Sciences,Key Laboratory of Biology and Genetic Improvement of Oil Crops,Ministry of Agriculture and Rural Affairs,Wuhan 430062,Hubei,China

Keywords:Brassica napus Silique length GWAS RNA-seq Phytohormones Cell number Cell size BnaA9.ARF18 BnaA9.CYP78A9

ABSTRACT Silique length influences seed yield in oilseed rape.It shows extensive variation in germplasm resources,and identifying the underlying genes and regulatory mechanisms would advance breeding for the trait.In the present study,a genome-wide association study (GWAS) using 331 core accessions planted in 10 environments revealed 13 loci associated with silique length on chromosomes A01,A04,A07,A09,and C03,explaining 6.2%–19.2% of phenotypic variance.Physiological analysis showed that silique length variation was attributable to differences in silique growth rate and/or duration before four weeks after flowering,with levels of endogenous phytohormones(auxin,ethylene,and GA24,GA12,and GA44)playing an important role.Cytological analysis showed that silique length variation was due mainly to differences in cell number followed by cell size.Transcriptomic analysis of two pools of silique walls with opposite length extremes revealed 3248 differentially expressed genes(DEGs).These DEGs were enriched in several pathways (such as cell wall,cell division,and hormone metabolism) associated with cell proliferation and expansion and silique development.Integrating GWAS,RNA-seq,and functional annotation results revealed 15 candidate genes for the major associated locus qSL.A09-3.Of these,BnaA9.ARF18 and BnaA9.CYP78A9 were validated by haplotype analysis followed by candidate gene association.Sequence variation in the coding region of BnaA9.ARF18 and expression of BnaA9.CYP78A9 in silique walls were strongly associated with silique length.Our results provide an explanation for the natural variation of silique length in oilseed rape germplasm and offer useful information for its improvement.

1.Introduction

Silique length is closely associated with seed yield in oilseed rape,being positively correlated with seed number per silique and seed weight [1–9].In oilseed rape germplasm resources,silique length shows large variation from several to more than 10 cm [3,4,8].However,the regulatory genes and mechanisms underlying this variation are poorly known.

From a physiological and developmental viewpoint,silique length is determined by growth rate and duration [10].However,the relative influence of growth rate and duration on its natural variation has not been studied.Fruit growth and development depend on endogenous phytohormones in the silique wall [11].However,the relationship between phytohormone levels and silique length is also ambiguous and invites investigation.At the cellular level,silique length is determined by proliferation and expansion of cells in the silique wall [11].The relative influence of cell number and size on silique length variation has not been investigated.

Genetically,silique length is a typical quantitative trait with continuous distribution,which is controlled by multiple genes and is readily affected by environmental conditions.More than 100 QTL for silique length have been identified in Brassica species by linkage and/or association analysis,and major QTL have been found on chromosome A09 [1–9,12–14].However,only two of these have been cloned and functionally characterized [10,15].

Reverse genetics studies have identified several genes that affect silique length in other Brassica species.Overexpression of two auxin-repressed protein genes BrARP1 and BrDRM1 in B.rapa reduced vegetative and reproductive growth,including the number and length of siliques,possibly via inhibition of either cell elongation or expansion [16].Overexpression of B.oleracea BoMF2 (a nuclear-localized AT-hook DNA-binding protein) in Arabidopsis led to shorter siliques and reduced pollen viability [17].Overexpression of the BcRISP1 gene in Arabidopsis may interrupt the mitochondrial electron transport chain by affecting the expression of mitochondrial respiration chain-associated genes,resulting in small siliques with reduced seed set [18].More than 100 genes that regulate/affect silique length have been identified in Arabidopsis and other plant species [11],of which most were cloned from mutants and unable to explain natural variation.

In the present study,we aimed to identify key genes and regulatory mechanisms underlying the natural variation of silique length in oilseed rape germplasm.By using an association population and the accessions with extreme long or short silique,the genetic,physiological,cytological,and molecular causes of silique length variation were investigated.

2.Materials and methods

2.1.Plant materials and field trials

An association population was constructed from 331 diverse oilseed rape accessions[19].All of them were grown in Nanchang,Jiangxi province in 2013 (N13) and 2014 (N14),in Wuhan,Hubei province from 2012 to 2017 (W12 to W17),and in Zhengzhou,Henan province in 2013 (Z13) and 2014 (Z14).The field planting followed a randomized complete block design with two replications.There were three rows 2 m long with 33 cm spacing in each plot,with 15 plants per row.Field management followed standard agriculture practice.At maturity,10 representative plants from the middle of each plot were harvested for silique length measurement following Li et al.[4].

2.2.Genome-wide association study

A genome-wide association study(GWAS)of silique length variation was performed using the association population.This population was previously genotyped using the Brassica 60 K Illumina R Infinium SNP array[19].The population structure was inferred with STRUCTURE 2.3.4 software[20].Three independent runs were performed with K(the putative number of clusters)varying from 1 to 10,and the length of burn-in period and number of Markov Chain Monte Carlo replications after burn-in were both set to 10,000.The most likely K value was determined from an ad hoc statistic ΔK based on the rate of change in LnP(D)between successive K values [21].Estimated membership coefficients for individual accessions (Q matrix;Table S1) were obtained from the corresponding runs.The relative kinship matrix (K matrix;Table S2) of the 331 accessions was calculated with TASSEL 5.2.41 software[22].

GWAS analysis employed the MLM (Q+K) model incorporated into TASSEL 5.2.41 software [22].Manhattan plots were constructed with R software [23].Linkage disequilibrium was estimated as the coefficient of correlation (r2) between allele’s frequencies of pairwise SNPs on each chromosome.The threshold for significantly associated SNP markers was set to P ＜4.08 10-5(P=1/24508,-log10(P)=4.39) as previously described [19].

2.3.Silique growth dynamics during development

To assess silique-length differences between long-and shortsilique accessions,silique length was recorded at 10 stages:0,3,6,9,12,15,18,21,24,and 28 days after flowering (DAF) and at maturity.To ensure the developmental consistency of the samples used for measurement,the flowers on the middle part of the main inflorescence and opening on the same date were tagged.At each stage,one flower/silique was taken from the tagged flowers of 10 representative plants.For each accession,10 flowers or siliques were used to measure ovary or silique length [26].

2.4.Quantification of phytohormones

To further investigate the physiological basis of silique length variation,the contents of endogenous phytohormones in the silique walls were measured at the key stage for the formation of silique length difference between the long-and short-silique accessions.At 15 DAF,the silique walls of long-and short-silique accessions were sampled and mixed equally to form two pools with three replications.The contents of six phytohormones including abscisic acid (ABA),Indole acetic acid (IAA),ethylene (ETH),Brassinosteroid (BR),cytokinin (CK) and gibberellin (GA) were measured on the phytohormone platform (http://www.genetics.ac.cn/jspt/zwjs/) of the Institute of Genetics and Developmental Biology,Chinese Academy of Sciences (Beijing,China).

2.5.Microscopic observation of cell number and size in ovary or silique wall

Cell number and size were measured in 12 long-and 11 shortsilique accessions at the initial stage of silique growth(the flowering date,0-DAF) and the stable stage of silique elongation,four weeks after flowering(4-WAF).For cytological analysis,four ovaries or siliques representing the mean length of each were sampled,and ovary or silique length,cell number,and mean length of cells in a long-axis section of the ovary or silique wall were recorded.The detailed experimental procedures for the section were described previously [9].

2.6.RNA sequencing and transcriptome analysis

To investigate how silique length differences form at the molecular level,RNA sequencing (RNA-seq) and transcriptome analysis were performed using the two contrasting pools.Silique walls at 15 DAF were sampled separately for the 12 long-and 11 shortsilique accessions chosen from the association population and then equally mixed to form two pools.RNA extraction,purification,and quantification were performed with the Plant Mini RNeasy Kit(Qiagen,Düsseldorf,Germany).RNA-seq libraries were prepared using the Illumina TrueSeq RNA Sample Preparation Kit (Illumina Inc.,San Diego,CA,USA).The final cDNA libraries were sequenced using the Illumina HiSeq platform,which generated raw reads.Low-quality,low-complexity,and repetitive raw reads were removed,and clean reads that passed quality control were used for subsequent analysis.All reads of each library were mapped to the reference genome of Darmor[24],and uniquely mapped reads were selected for transcript quantification.False discovery rate(FDR) 0.05 and P-value 0.005 were used as a threshold to assign differentially expressed genes (DEGs).The raw RNA-seq data have been deposited in NCBI Sequence Read Archive (SRA,http://www.ncbi.nlm.nih.gov/sra/) as accession GSE176120.

2.7.Validation of differentially expressed genes

Using the same RNA samples for sequencing as template,reverse transcription was performed using the PrimeScriptTMII 1st Strand cDNA Synthesis Kit (Takara,Osaka,Japan),according to the manufacturer’s protocol [25].Quantitative real-time PCR(qRT-PCR) was performed using Bio-Rad IQ5 with SYBR Green detection (https://www.bio-rad.com/).Ten DEGs were randomly selected and validated by qRT-PCR using gene-specific primers designed by Primer Premier 5.0 software (http://www.premierbiosoft.com/).The β-actin gene in oilseed rape was used as an internal control to normalize the expression levels.A cycling temperature of 57°C together with a single peak on the melting curve was used as the standard to confirm the specificity of designed primer pairs,and each sample was assessed in triplicate [9].

2.8.Editing the orthologs of candidate genes in Arabidopsis thaliana

The coding sequences (CDS) of candidate genes were subjected to BLAST analysis against the Arabidopsis genome (https://www.arabidopsis.org/Blast/index.jsp) to identify orthologs.The orthologs in Arabidopsis (Col-0) were targeted for gene editing to verify their function.Search software (http://crispr.hzau.edu.cn/CRISPR2/)was used to identify small guide RNA (sgRNA) sequences within the CDS of the target genes (near the start codon) and determine the optimal sgRNAs according to the on-target efficiency scores,and possible off-target sites [27].

Escherichia coli strain DH5α and Agrobacterium tumefaciens strain GV3101 were used for vector construction.AtU6-26-sgRNA-SK and pCAMBIA1300 were the intermediate and final vectors,respectively,for sgRNA cassette construction and plant transformation.Gene knockout was performed with the CRISPR/Cas9 system following Yan et al.[28]

2.9.Statistical analysis

Descriptive statistics,analysis of variance (ANOVA),and multiple comparisons were performed using SAS version 8.0 software(https://www.sas.com/en_us/home.html).Broad-sense heritability was calculated as:h2=Whereis the variance of genotype,is the variance of the genotype-byenvironment interaction,is the error variance,n is the number of environments,and r is the number of replications.

3.Results

3.1.Genetic basis of silique length variation in oilseed rape germplasm

The silique lengths of the 331 accessions in 10 environments varied from 29.0 mm to 108.4 mm(3.73 fold),with the coefficient of variation ranging from 0.13 to 0.20 (Table S3).As expected,the silique lengths of the 331 accessions showed normal or nearnormal distributions in all 10 environments (Fig.1A).Analysis of variance showed that genotype,environment,and their interaction exerted significant effects on silique length (Table S4).The broadsense heritability of silique length was as high as 0.96,suggesting its high stability suitable for genome-wide association study.

Of 24,508 SNPs,56 were significantly associated with silique length(Fig.1B;Table S5).Of these,1,4,1,1,9,6,6,8,4,and 16 were from N13,N14,W12,W13,W14,W15,W16,W17,Z13,and Z14,respectively.Removal of closely related SNPs within 500 kb left 25 nonredundant SNPs (Table S6).These were further combined into 13 associated loci,which were distributed on chromosomes A01,A04,A07,A09 and C03,explaining 6.6%–19.2% of the phenotypic variance (Table 1).Two associated loci on chromosome A09 (qSL.A09-2 and qSL.A09-3) could be repeatedly identified in multiple environments.qSL.A09-3 was detected in nine of the ten environments,explaining 11.2%–19.2% of the phenotypic variance,and qSL.A09-2 was detected in five environments,explaining 6.7%–17.7% of the phenotypic variance.As expected,the two peak SNPs(Bn-A09-p29953651 and Bn-A09-p28996871)were strongly associated with silique length variation in the association panel(Fig.S1A,B).The remaining 12 SNPs were identified in only one environment and explained 6.6%–9.1%of the phenotypic variance(Table 1).

Comparative analysis based on physical position on the genomic map of Darmor_V4.1 (https://www.ncbi.nlm.nih.gov/assembly/GCA_000751015.1) showed that three SNPs identified in this study overlapped reported QTL (Table S7),whereas the remaining 10 loci appear to be novel (Table 1).

To further dissect the regulatory mechanisms underlying the natural variation of silique length in this association population,representative accessions with extremely long or short siliques were chosen for a systematic comparative study in physiological,cellular and molecular levels.

3.2.Physiological basis of silique length variation

3.2.1.Silique length variation is caused by differences in growth rate and/or duration

At the flowering date (0-DAF),the ovary lengths of long-and short-silique accessions ranged from 4.2 to 10.3 mm and 3.7 to 8.2 mm,respectively,and the mean of the long (7.9 mm) was greater (P=1.6E-02) than that of the short (6.1 mm),with a 1.3-fold difference (Fig.2).However,at 4-WAF,the silique lengths of the long-and short-silique accessions increased from 67.2 to 90.6 mm and 33.6 to 45.5 mm,respectively,and the mean of the long (80.7 mm) was greater (P=6.2E-13) than that of the short(40.0 mm),with a 2.0-fold difference.Thus,the increase in mean length for the long-silique accessions (10.2-fold) was greater than that (6.6-fold) for the short-silique accessions.The silique lengths of long-and short-silique accessions increased rapidly from 12 to 15 DAF and 9 to 12 DAF,respectively,a difference of three days.The stable/fix stage of silique length for the short-silique accessions was earlier than that for long-silique accessions,suggesting that the growth duration of long-silique accessions was greater than that for short-silique accessions.The long-silique accessions generally grew faster than the short-silique accessions.Using the silique growth characteristics described above,the long-silique accessions could be divided into fast and sustained growing types.Thus,the final difference in length between the extremely longand short-silique accessions was caused by the difference in its growth rate and/or duration.

3.2.2.Silique length variation was associated with endogenous phytohormone levels

Levels of six endogenous phytohormones (ABA,IAA,ETH,BRs,CK,and GA) were measured in the silique walls of long-and short-silique accessions (Table 2).The contents of IAA and ETH in the silique walls of long-silique accessions (161 and 468 pg mg-1)were higher than in those of short-silique accessions (150 and 329 pg mg-1).In contrast,the contents of ABA,all detected BRs and CKs,and three bioactive GAs(GA1,GA3,and GA4)were lower in silique walls of long-silique than in those of short-silique accessions.One type of BR(TE),two types of CK(DHZ and iP9G),and GAs(GA7 and GA29) were not detected in silique walls of either longor short-silique accessions,while one GA (GA3) was undetected only in long-silique accessions.

3.3.Cytological basis of silique length variation

At the flowering date (0-DAF),the ovary lengths of 12 longsilique and 11 short-silique accessions used for sectioning ranged from 5.5 to 8.9 mm (1.62-fold) and 4.1 to 7.3 mm (1.78-fold),respectively (Fig.S2A),and the mean (7.2 mm) of the former was greater (1.2-fold;P=1.53E-02) than that of the latter (6.0 mm).The numbers of cells in the long-axis section of the silique walls of the 12 long-and 11 short-silique accessions ranged from 544.5 to 940.5 (1.73-fold) and 365.0 to 704.0 (1.93-fold),respectively (Fig.3A),and the mean of the former (710.0) was greater(P=7.97E-03) than that of the latter (566.4).However,the calculated mean lengths of cells in the long-axis section of the silique wall varied from 8.0 to 11.7 μm (1.46-fold) and 8.5 to 11.9 μm(1.40-fold),respectively,for the 12 long-and 11 short-silique accessions and their means (10.3 and 10.7 μm) were not different(Fig.S2B).

Fig.1.GWAS of silique length in 10 environments.(A)Frequency distribution of silique length of the association population planted in 10 environments.The horizontal axis shows the silique length grouped by 5-mm increments from 25 to 110 mm,and the vertical axis shows the number of accessions in each group.(B)Manhattan plot for GWAS of silique length using 331 accessions planted in 10 environments.The horizontal axis shows the lengths of the 19 chromosomes,and the vertical axis shows the–log10(P)value for each SNP marker;the dotted line shows the significance threshold.The colors correspond to environments (Nanchang 2013–2014,Wuhan 2012–2017,and Zhengzhou 2013–2014).

At 4 WAF,the silique length of 12 long-and 11 short-silique accessions ranged from 65.0 to 91.6 mm and 34.1 to 43.9 mm,respectively(Fig.S3A),and the mean of the former(78.4 mm)was greater(2.05-fold;P=1.09E-12)than that of the latter(38.2 mm).The number of cells in the long-axis section of the silique wall ranged from 2307 to 2834 and 1084 to 1599 respectively for the long-and short-silique accessions,and the mean of the former (2527) was greater (1.87-fold,P=4.65E-13) than that of the latter (1350)(Fig.3B).The mean lengths of cells in the long-axis section of the silique wall varied from 28.2 to 35.9 μm and 24.9 to 37.0 μm respectively,for the long-and short-silique accessions.Although the mean of the former (31.0 μm) was higher than that of the latter(28.3 μm),the difference was not statistically significant(Fig.S3B).

3.4.Comparison of transcriptomes of long-and short-silique accessions

3.4.1.Transcriptome sequencing and mapping

After filtering,respectively 209,313,464 and 216,867,458 clean reads were obtained from RNA sequencing of three repeats of long-and short-silique walls.The mean Q30 percentages for the three repeats of long-and short-silique pools were 90.9%.The guanine and cytosine(GC)contents for the three repeats were respectively 46.3% and 46.5%.Most of the clean reads were successfully mapped to the reference genome,with mean mapped rates of 90.7% and 90.9% for the long-and short-silique pools (Table S8).The unique-and multiple-mapped rates of the six samples varied from 60.9% to 68.1% and from 22.5% to 30.1%,with means of 66.2% and 24.6%.Of the 101,040 genes detected in the two pools,3248 differentially expressed genes(DEGs)were identified,including 1116 up-regulated and 2132 down-regulated genes (Fig.S4A,B).Although the fold changes of the 10 selected DEGs in RNA-seq(Table S9) were generally higher than those in qRT-PCR experiments,the overall trend of up-and down-regulation was consistent between the two methods (Fig.S5),supporting the reliability of the DEGs.

3.4.2.DEGs in the silique wall were strongly associated with silique development

The 3248 DEGs (Table S10) were classified using SuperViewer(http://bar.utoronto.ca/ntools/cgi-bin/ntools_classification_superviewer.cgi),showing the total number and frequency as well as theP-value for each class (Fig.S6A,B).Of the 21 enriched (P ＜0.05)classes,11 reached the highly significant level (P ＜0.001),including DNA,hormone metabolism,RNA,misc,not assigned,cell wall,transport,development,secondary metabolism and protein,reflecting the different developmental and metabolic statuses of the two pools.

Table 1 Details of 13 peak SNPs associated with silique length.

Fig.2.Silique growth dynamics of representative long-and short-silique accessions.The vertical axis shows the silique length in millimeters,and the horizontal axis shows different stages of silique growth.Colors correspond to accessions.

Of the 284 DEGs grouped into RNA class (Table S8),the largest class,251(88.4%)were involved in regulation of transcription.This result demonstrated the importance of transcription factors (TFs)in regulating the expression of genes involved in silique development.These TFs were from nearly 20 families,including zinc finger(41),MYB (26),AP2/EREBP (22),bHLH (21),WRKY (12),NAC (11),Homeobox (10),B3-domain (9),bZIP (7),Trihelix (6),MADS (6),AUX/IAA(6),HSF(4)and TCP(4)etc.Of the 239 DEGs grouped into protein class,only five were involved in amino acid activation and protein folding;18 were involved in protein targeting.In contrast,the other 39,57 and 112 DEGs belonged to the sub-classes protein synthesis (most as ribosomal protein),posttranslational modification,and degradation (most as E3 ubiquitin ligase and protease).The 161 DEGs grouped into the misc class,except for 10 lipid transfer proteins (LTP),belonged to various enzymes,including cytochrome P450 (30),peroxidases (17),UDP glucosyl and glucoronyl transferases (16),short-chain dehydrogenase/reductase (7),GDSL-motif lipase(9),gluco-,galacto-and mannosidases(9),nitrilases(10),glutathione S transferases(11)and oxidases-copper(7),myrosinases-lectin-jacalin(7)and acid and other phosphatases(6).A total of 110 DEGs were grouped into the transport class,associated with sugars (18),ABC transporters and multidrug resistance systems(16),metabolite transporters(11),peptides and oligopeptides (10),p-and v-ATPases (7),and nucleotides (6).A set of 83 DEGs were involved in the metabolic pathways of several phytohormones,including AUX (19),GA (17),ETH (16),ABA (11),JA(7),BR (6),CK (5) and SA (2),demonstrating the complex role and interplay of phytohormones in silique development.A total of 73 DEGs were grouped into cell wall pathways and could be divided into several subclasses,including precursor synthesis(6),cellulose(9) and hemicellulose (3) synthesis,cell wall proteins (11) mainly as AGPs,degradation related mannan-xylose-arabinose-fucose,pectate lyases and polygalacturonases (12),modification (16) and pectin esterases (13).A total of 57 DEGs belonged to secondary metabolism pathways,including flavonoids (14),glucosinolates(12),simple phenols (11),phenylpropanoids (10),isoprenoids (6),and alkaloid-like (3).Of the 31 DEGs grouped into the DNA class,respectively 21 and 4 were involved in synthesis/chromatin structure and repair,which provide a structural basis for cell expansion and proliferation.

Table 2 Endogenous phytohormone contents in silique walls of long-and short-silique accessions.

To establish a link between these DEGs and silique development,they were compared with reported fruit size genes,mainly from Arabidopsis [11].More than 10 DEGs were similar to known fruit size genes,including CYP78A9 (BnaA04g00510D,BnaA09g55530D,BnaC08g31760D),CYP72C1 (BnaCnng42960D),BT5(BnaA03g54130D and BnaC07g46630D),AGL63 (BnaA09g25390D),AGP19 (BnaA07g27370D),AHP4 (BnaC03g39610D),CBSX4(BnaCnng33730D),GID1B (BnaA09g55590D),NMNAT(BnaA03g11210D),and RMI1 (BnaA09g06650D).More than half of the DEGs were involved in genetic and signaling pathways of fruit size regulation[11],including the ubiquitin–proteasome pathway,phytohormones,transcription factors,receptor kinases,G-protein signaling,arabinogalactan,RNA-binding proteins,cell organization,and cell wall proteins (Table S10).

3.5.Identification and verification of candidate genes for the major loci

3.5.1.Identification of candidate genes for qSL.A09-3 by integration of GWAS and RNA-seq with functional annotation

Although 13 associated loci were identified,only qSL.A09-3 was found in almost all the environments.It also explained the largest variance from 10.8%to 19.2%,and was accordingly chosen as a target for further research.By the integrative analysis of GWAS and RNA-seq with functional annotation,a total of 15 candidate genes were identified for qSL.A09-3(Table 3),including two known genes for silique length [10,15],BnaA09g55580D (BnaA9.ARF18) and BnaA09g55530D (BnaA9.CYP78A9).BnaA09g55530D showed the most significant difference(P=3.13E-248)among the 3248 DEGs,whereas BnaA09g55580D was absent from the list of DEGs(Table S10).Of the 15 candidate genes,BnaA9.ARF18 and BnaA9.CYP78A9 were validated only by candidate gene association analysis,whereas the other 13 were preliminarily verified by editing their orthologs in Arabidopsis.

Table 3 List of 15 candidate genes identified for the major silique length-associated locus qSL.A09-3.

3.5.2.Validation of BnaA9.ARF18 and BnaA9.CYP78A9 by candidate gene association analysis

3.5.2.1.Haplotype analysis of BnaA9.ARF18 and BnaA9.CYP78A9.Because in the transcriptomic analysis of silique walls,the expression of BnaA9.ARF18 showed no significant difference between the long-and short-silique accessions,sequence and haplotype analysis of BnaA9.ARF18 was limited to the coding sequence.The coding sequences of the BnaA9.ARF18 gene in the 48 representative oilseed rape accessions could be classified into nine haplotypes,based on nucleotide and amino acid sequence variation (Figs.S7A,S8).The frequencies of the nine haplotypes showed large differences,varying from 1 (2.3%) to 18 (40.9%) accessions (Fig.S7B).In particular,haplotypes seven and five were present in 18 and 17 accessions,which accounted for the highest proportions respectively 40.9%and 38.6%.Among the 16 accessions with long siliques,the BnaA9.ARF18 sequence of 12 accessions belonged to haplotype seven,with a proportion 75.0%that was much higher than the corresponding one(40.9%)in all the investigated accessions(Fig.S7C).In addition,of the 16 accessions with short silique,the BnaA9.ARF18 sequence of 13 accessions belonged to haplotype five,showing a proportion(81.3%)that was also much higher than the corresponding one (38.6%) in all the investigated accessions (Fig.S7D).These results strongly suggested that haplotypes seven and five of BnaA9.ARF18 were associated with respectively the long and short silique phenotypes.

The full-length genomic DNA of BnaA9.CYP78A9 (including upstream promoter,transcribed region,and 3′-downstream) in the 48 representative oilseed rape accessions were also cloned and sequenced.Only the upstream promoter(rather than the transcribed region and 3′-downstream) showed sequence variation among these accessions,which could be classified into five haplotypes (Fig.S9A).Among all the 48 investigated accessions(Fig.S9B),haplotype five was present in 21 accessions (43.8%),haplotype one in 15 accessions (31.3%),haplotype two in seven accessions (14.6%),haplotype three in four accessions (8.3%),and haplotype four in only one accession (2.1%).Among the 16 accessions with long silique (Fig.S9C),the BnaA9.CYP78A9 sequence in 15 accessions belonged to haplotype five,showing a proportion of 93.8% that was much higher than the corresponding one(43.8%) in all the investigated accessions.For the 16 accessions with short silique (Fig.S9D),the BnaA9.CYP78A9 sequence of 11 accessions belonged to haplotype one,with a proportion (68.8%)that was also higher than the corresponding one (31.3%) in all the investigated accessions.These results strongly suggested that haplotypes five and one of BnaA9.CYP78A9 were associated with the long-and short-silique phenotypes,respectively.

Fig.3.Cytological basis for the natural variation of silique length.Number of cells in the long-axis section of the inner epidermis for(A)0-DAF ovary and(B)4-WAF silique wall.The horizontal axis shows the 12 long-silique and 11 short-silique accessions,and the vertical axis shows cell numbers.The heights of columns and bars on columns represent mean and standard deviation,respectively.Lowercase letters following values on bars identify values that are not significantly different in multiple comparisons.

3.5.2.2.Association of BnaA9.ARF18 and BnaA9.CYP78A9 expression level with silique length.The relative expression levels of BnaA9.ARF18 in the 12-DAF silique walls of the 48 representative oilseed rape accessions showed moderate variation from 0.26 to 1.46.As expected,the mean expression levels of BnaA9.ARF18 in the silique walls of the long-and short-silique accessions were 0.64 and 0.87,respectively,which were not significantly different(Fig.S10A).The correlation(r=–0.36)between the expression level of BnaA9.ARF18 in the silique walls and the silique lengths of the 48 investigated accessions was also small and negative(Fig.S10B),which accorded well with its negative regulation of silique development [15].

In contrast,for BnaA9.CYP78A9 (Fig.S10C),its relative expression levels in the 12-DAF silique walls of the 48 investigated accessions showed significant(1339-fold)variation from 0.11 to 147.32.The expression levels of BnaA9.CYP78A9 were generally higher in the silique walls of the long-silique than in those of the shortsilique accessions,and the mean (83.29) of the former was 20.8-fold that of the latter (4.01).As expected,the expression level of BnaA9.CYP78A9 in the silique walls was also highly positively correlated(r=0.68)with the silique length of the investigated accessions (Fig.S10D),which also accorded with its positive regulation of silique development [10].

3.5.2.3.Association of BnaA9.ARF18 and BnaA9.CYP78A9 sequence variation with silique length.To further confirm that BnaA9.ARF18 and BnaA9.CYP78A9 were causal genes for the natural variation in silique length,sequence variations between haplotypes seven and five of BnaA9.ARF18 and between haplotype five and one of BnaA9.CYP78A9 were used to develop SNP markers.One SNP marker in the coding region of BnaA9.ARF18 and another in the promoter region of BnaA9.CYP78A9 were developed and used for the genotypic analysis of 331 accessions in the association panel.

For the BnaA9.ARF18 gene,the two variant bases A and G,accounted for respectively 84.7% and 15.3% of all the 331 accessions in the association panel.The G allele’s mean silique length(73.46 mm) was much greater than that (58.3 mm) of the A allele(Fig.S11A).

For the BnaA9.CYP78A9 gene,the two variant bases A and G accounted for respectively 57.0% and 43.0% of all the 331 accessions in the association panel.However,the mean silique length(55.61 mm) of the G allele was slightly higher than that(61.14 mm) of the A allele (Fig.S11B),suggesting that this SNP might not be the functional variation in this gene.

3.5.3.Validation of other candidate genes by editing of their orthologs in Arabidopsis thaliana

The AT3G61460/BRH1-edited lines showed shortened siliques and decreased numbers of seeds (Fig.4).The mean silique length of four edited lines(M1–M4)was 5.13 cm,65.4%of that of the wild type(7.85 cm in W1–W2).The mean seed number per silique in the four edited lines was 21.5,which was also much lower than that(34.0)of the wild type.These results confirmed that the nucleotidesequence of the target site was edited efficiently (Fig.S12).The BRH1 gene encodes a novel ring finger protein with ubiquitinprotein ligase activity in the brassinosteroid(BR)mediated signaling pathway.Overexpression of BRH1 in Arabidopsis led to rounder rosette leaves [29].

4.Discussion

4.1.Novel loci identified for silique length

Among the 13 association loci identified in the present study,three were consistent with previously published silique length QTL,whereas the other ten appear to be novel (Table 1).Two repeatedly identified associated loci (qSL.A9-2 and qSL.A9-3) on chromosome A09 showed a large effect and were consistent with the previous published major QTL for silique length [2,3,4,5,8,9].As expected,the two peak SNPs underlying major associated loci,Bn-A09-p28996871 and Bn-A09-p29953651,were also strongly associated with silique length variation in the association panel(Fig.S1A,B),and should be used in marker-assisted selection procedures.To date,only two major silique length QTL on chromosome A09 have been cloned [10,15],and both causal genes(BnaA9.ARF18 and BnaA9.CYP78A9) are very close to the peak SNP of qSL.A9-3 (Bn-A09-p29953651).Another repeatable major association locus,qSL.A9-2,remains to be cloned.

4.2.Physiological basis of silique length variation in oilseed rape

Observation of silique growth dynamics showed that final silique length difference was determined mainly by the silique growth process,followed by the initial ovary length difference(Fig.2).This conclusion is consistent with a previous finding [30]that fertilization triggers fruit development in plants.The silique length difference between long-and short-silique accessions was attributable to both growth rate and duration.Some long-silique phenotypes (such as in accessions 3S1145 and 3S1129) were caused mainly by fast growth,whereas others (such as in accessions 3S1202 and 3S1214)were due mainly to long duration.It should thus be possible to create super-long-silique accessions by combining the two characters of fast growth and long duration.In oilseed rape,seed filling follows the completion of silique morphogenesis,usually about 25 days after flowering[31].For this reason,fast silique growth should be favorable (especially for earlymaturing cultivars) because it leads to early seed filling and oil accumulation.The finding that the highest growth rates occurred at different stages in different cultivars may be attributed to genetic factors and phytohormone levels and their interactions with environmental conditions.

IAA and ETH were more abundant in the silique walls of longsilique accessions.In contrast,ABA,BR,CK,and three bioactive gibberellins(GA1,GA3,and GA4)were more abundant in short-silique accessions (Table 2).These results suggest that differences in the six phytohormones contribute to silique length variation in oilseed rape.Consistent with this suggestion,of the 83 DEGs involved in phytohormone pathways,35 were involved in the synthesis or degradation of these phytohormones,thereby presumably contributing to the phytohormone abundance difference between the long-and short-silique walls.All three identified genes for silique length involved the phytohormone pathway:BnaA9.CYP78A9[10]and BnaA9.ARF18[15]are involved in respectively the synthesis and signal transduction of AUX,while BnaA9.BRH1 is involved in the signal transduction of BR.These results strongly suggest the importance of phytohormones in regulating silique length.

The further cytological analysis showed that the silique length difference between the long-and short-silique accessions was due mainly to the difference in cell number (Fig.3),followed by cell size(Fig.S3).However,both cloned major QTL genes in oilseed rape affect cell size rather than cell number [10,15].The disagreement may be due to the difference in the tissues sampled for observation:in the present study,the endocarp was used for sectioning,whereas in the above two studies,the pericarp was directly used for observation by scanning electron microscope.In a previous study[32]in Arabidopsis,different trends in changes in cell number or size were observed in different layers of the silique walls of global and ga1 mutants and a wild type.

4.3.Identification of causal genes using an integrated strategy

A review of cloned plant QTL stated that they could result from changes in protein function and/or expression level [33].An analysis integrating GWAS,RNA-seq,and functional annotation identified only 15 candidate genes for the major locus qSL.A9-3(Table 3).Two strategies were used to validate these candidate genes:for those known to regulate silique length,association analysis was used [34];whereas for other candidates,the gold standard of a transgenic experiment was used.Using 48 representative oilseed rape accessions selected from a previously developed association population in our laboratory [9],the sequence variations of BnaA9.ARF18 and BnaA9.CYP78A9 [10,15] were revealed first at the population level (Figs.S7A,S9A).

The SNP markers developed from the coding sequence of BnaA9.ARF18 were strongly associated with silique length in the association panel (Fig.S11A).In contrast,the relative expression levels of BnaA9.ARF18 in silique walls did not differ between the longsilique and short-silique accessions and were not correlated with their silique lengths(Fig.S10B).These results strongly suggest that the functional variation of BnaA9.ARF18 is located in the coding region,a conclusion in accord with the previous functional characterization of this gene [15].However,no sequence variation was found in the transcribed region of BnaA9.CYP78A9,and an SNP marker developed from its promoter region was weakly associated with silique length in the association panel(Fig.S11B).The relative expression levels of BnaA9.CYP78A9 in the silique walls of representative oilseed rape accessions were highly positively correlated(r=0.68)with silique length(Fig.S10D),a finding also consistent with its positive regulatory role in silique development[10].These results suggest that the functional variation of BnaA9.CYP78A9 is associated with its expression abundance,a conclusion in accord with the recent finding that a CACTA-like transposable element inserted in the upstream regulatory region of BnaA9.CYP78A9 acted as an enhancer to increase its expression [10].

Fig.4.Editing of brh1 gene in Arabidopsis.(A)Representative plants and(B)siliques of wild type and brh1 edited lines.(C)Silique length and(D)seed number per silique of representative wild type (W1–W2) and brh1-edited lines (M1–M4).

We identified a new gene affecting silique length,BnaA09g39200D,by editing its ortholog BRH1 in Arabidopsis.Several other genes in BR pathway are known to regulate silique length [11],such as BRI1 and SHK1.In addition to reducing silique length (Fig.4C),knockout of this gene also reduced seed number per silique (Fig.4D).This finding was in accord with our previous finding [11] that a change of seed number per fruit was usually accompanied by a change of fruit length in the same direction.Overexpression of BRH1 in Arabidopsis led to a change in leaf shape[29].These results suggest that BRH1 plays multiple roles in plant growth and development.Future studies should verify the function of BnaA09g39200D in Brassica,mainly to clarify whether its overexpression could increase silique size and seed yield of oilseed rape.

4.4.A proposed model for silique length variation

We propose a silique-length regulatory model:variation in the sequence and/or expression of causal genes (e.g.,BnaA9.CYP78A9)underlying associated loci for silique length →change in endogenous phytohormone level in the silique wall (e.g.,AUX) →perturbation of signal transduction pathways(e.g.,receptor kinase signaling)→activation or repression of the expression of downstream responsive genes(e.g.,IAA28)→promotion or repression of the division and/or growth of cells in the silique wall →change in the number and/or size of cells in the silique wall →variation in silique length.

In summary,our study identified key genes and regulatory mechanisms underlying the natural variation of silique length in oilseed rape,providing a theoretical basis and technical support for its genetic improvement.

CRediT authorship contribution statement

Quaid Hussain:Data curation,Investigation,Visualization,Writing– original draft.Jiepeng Zhan:Investigation.Huabing Liang:Validation.Xinfa Wang:Resources.Guihua Liu:Resources.Jiaqin Shi:Conceptualization,Project administration,Investigation,Methodology,Writing– review &editing.Hanzhong Wang:Formal analysis,Funding acquisition,Supervision.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

We thank the anonymous reviewers for improving this manuscript.We also thank Prof.Xue-Rong Zhou (Senior Principal Research Scientist in CSIRO Agriculture &Food) who has carefully read the whole manuscript and given us many constructive comments and suggestions.This research was supported by the Technical Innovation Project of Hubei Province (2018ABA087),National Natural Science Foundation of China (31771840),the China Agriculture Research System of MOF and MARA (CARS-13),the Agricultural Science and Technology Innovation Program(CAAS-ASTIP-2013-OCRI),and Fundamental Research Funds for Central Non-Profit Institute of Crop Sciences CAAS (Y2020YJ09).

Appendix A.Supplementary data

Supplementary data for this article can be found online at https://doi.org/10.1016/j.cj.2021.08.010.