The DH population consisting of 201 lines was originated from two Spring and semi-Winter rapeseed varieties 3S1181 and 3S1235, which were chosen from an association population constructed in our lab (Qadir et al., 2022). Field experiments were conducted in the typical spring and semi-Winter rapeseed planting areas across three years, including Ping’an in Qinghai province (32.27° N, 114.52° E) in 2020 (coded as 20QH), Yangluo in Wuhan city (29.58° N, 113.53° E) in 2021 and 2022 (coded as 21YL and 22YL), and Huanggang city in Hubei province (30.44° N, 114.87° E) in 2022 (coded as 22HG). The field experiments were conducted according to the procedure described previously (Qadir et al., 2022).

For investigating ONPO, three buds on the main inflorescence of five plants were randomly sampled at the beginning of flowering, which resulted in 15 buds for each replication. The formalin solution was used to fix and store the samples buds (Yang et al., 2017). Petals and stamens were removed by tweezers and the ovary was taken out, according to the previously published experimental procedures (Qadir et al., 2022).

In order to examine whether the variation of ONPO is related to the levels of endogenous phytohormones, auxin (AUX), Brassinolide (BR), cytokinin (CK) and gibberellin (GA) in the ovaries (at the initial stage of ovule development) of two pools of high and low ONPO lines were measured. The buds between 0.5 to 1 mm length were selected for dissection by tweezers under microscope. The frozen ovaries from high or low ONPO lines in each biological repeat were equally mixed to produce the corresponding pools. The abundance of phytohormones was detected with phytohormonal platform at Metware Biotechnology Company (Wuhan, China).

The two parents and 201 DH lines were genotyped using the 50K SNP array comprising of 45,707 SNP markers, which resulted in a total of 11348 polymorphic markers. Moreover, after removing the markers with heterozygosity >5% or missing rate >20%, a total of 7176 high-quality SNP markers were obtained for subsequent analysis. The genetic linkage analysis was performed using JoinMap 4.0 software (Van Ooijen, 2006) with the default parameters: goodness of fit ≤5.0, LOD scores >2.0, recombination frequency < 0.4.

Linkage mapping of QTL was performed using CIM (composite interval mapping) method in WinQTLCart2.5 software with the default setting (Zeng, 1994). The LOD (logarithm of the odds) threshold for defining a significant QTL was determined by permutation test (Churchill and Doerge, 1994) with 1000 repetitions at P < 0.05. The LOD value, additive effect (A), confidence interval (CI), and phenotypic variance explained (PVE) of each QTL were estimated.

The dissected ovaries (at the initial stage of ovule development) from high or low ONPO lines were equally mixed to produce the corresponding pools with three biological repeats (M1-M3 and L1-L3). The total RNA of each sample was isolated using the RNeasy Plant Mini Kit (Qiagen Inc., Germany). The cDNA library was constructed and sequenced using DNBSEQ-T7 platform of BGI Genomics Co., Ltd (Shenzhen, China).

The low-quality raw reads were removed by filtering with SOAPnuke v1.4.0 software (Chen et al., 2018), and the remaining clean reads were used for subsequent analysis. The clean reads were mapped to Darmor-bzh reference genome using HISAT v2.1.0 software (Kim et al., 2015). The transcripts of each sample were re-constructed and integrated by StringTie and Cuffmerge software respectively (Shumate et al., 2022). The gene expression abundance in each sample was estimated using RSEM v1.12.8 software (Li and Dewey, 2011). Both the P-value 0.005 and Q-value 0.05 were simultaneously used to define the deferentially expressed genes (DEGs). The clean reads of the six samples were uploaded to NCBI database with accession number PRJNA1150986.

The classification and enrichment analysis of Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were performed using the online bioinformatics analysis platform of BGI Genomics Co., Ltd (https://biosys.bgi.com). The functional classification of DEGs were conducted by online tool Super-Viewer (http://bar.utoronto.ca/ntools/cgi-bin/ntools_classification_superviewer.cgi).

To validate the reliability of DEGs identified by RNA sequencing, qRT-PCR was performed using the same RNA for sequencing. The M-MLV reverse transcriptase (Promega) was used to synthesize cDNA containing total RNA (4μg) and oligo (dT) primers, as specified by the company instructions. The qRT-PCR experiment was performed using SYBR Green Real-time PCR Master Mix (Toyobo, Osaka, Japan) and a Bio-Rad CFX96 Real-Time Detection System. The relative expression was calculated using the 2−ΔΔCT method, with B. napus ACTIN2 as an internal control (Qadir et al., 2022).

The phenotypic data of the ONPO was statistically analyzed using Microsoft Excel, including mean, standard deviation, and coefficients of variation etc. The phenotypic correlation and variance were respectively analyzed using the CORR and ANOVA procedure in SAS software. The broad-sense heritability was estimated using the following formula: h ² = σ_g ²/(σ_g ² + σ_ge ^2/n + σ_e ²/nr), where σ_g ², σ_ge ² and σ_e ² respectively represented the variances caused by genotype, genotype×environment and error, while n and r respectively represented the environment and replicate number.

The two parents and 201 DH lines were planted in four environments including 20QH, 21YL, 22YL and 22HG. Of these, 20QH belongs to Spring type environment and seeds were sown at the middle of May whereas 21YL, 22YL and 22HG belong to semi-Winter type environment and seeds were sown at the end of September. The ONPO of parent 3S1235 in four investigated environments ranged from 32.7 to 36.8 in average, and was significantly more than that of another parent 3S1181 from 27.2 to 29.3 in average ( Table 1 ). In each of the four investigated environments, the ONPO range of DH population exceeded that of two parents, indicating that favorable alleles should be present in both parents. Interestingly, the ONPO mean (35.3) of DH population in 20QH was much higher than those in other three environments, suggesting that oilseed rape tended to differentiate more ovules in Spring than semi-Winter type growing areas. As expected, the ONPO of DH population showed normal distribution in all four investigated environments ( Figure 1 ), exhibiting quantitative genetic characteristics.

Generally, the ONPO of DH population in three semi-Winter rapeseed environments (21YL, 22YL and 22HG) were highly correlated (r=0.87, 0.85 and 0.85), but they exhibited low correlation with Spring rapeseed environment 20QH (r=0.39, 0.28 and 0.38) ( Supplemenmtary Table S1 ), which suggested that ONPO might vary across the different environments. The results of variance analysis demonstrated that genotype (G), environment (E), and their interaction (G × E) were all significant on ONPO ( Supplementary Table S2 ). The calculated broad-sense heritability for ONPO was as middle as 0.52, likely due to the large environmental difference in 20QH.

The remaining high-quality SNP markers were subjected to linkage grouping and resulted in a skeleton genetic map of 2111 SNP markers within 19 linkage groups ( Figure 2 ; Table 2 ). The number of SNP markers in 19 linkage groups varied from 68 (A07) to 171 (C04), with an average of 111. The genetic length of 19 linkage groups varied from 63.6 (A02) to 118.1 (C3) cM, with a total length of 1715.7 cM. The marker density of 19 linkage groups varied from 0.76 (C05) to 1.95 (A09), with a mean of 1.23 per cM. The physical length of 19 linkage groups varied from 15.8 (A10) to 54.4 (C3) Mb, with a total length of 587.9 Mb. The covered ratios of 19 linkage groups varied from 63.6% (C06) to 100.0% (C04), with a mean of 91.5%. It was obvious that the recombination frequencies of the 10 linkage groups in A sub-genome (2.63 to 5.36, mean=3.94) were mostly higher than the 9 linkage groups in C sub-genome (1.79 to 3.58, mean=2.35), with a mean of 2.92 per Mb for the whole genetic map.

A total of 10 consensus QTLs of ONPO were identified in the DH population across four investigated environments ( Figure 3 ; Table 3 ). Of these, 5, 3, 3 and 2 QTLs were identified from 20QH, 21YL, 22YL and 22HG environments, respectively. These QTLs were distributed in A01, A03, A06, A07, A09, A10, C02, and C05 linkage groups, accounting for 7.0-15.9% of the phenotypic variance. The additive effects of these QTLs could be either positive or negative, ranging from -1.16 to 1.32, which was highly accordant with the over-parent distribution of ONPO in DH population.

It should be noted that two QTLs in A07 and C02 linkage groups (qONPO.A07-1 and qONPO.C02-1) were repeatedly detected in different environments and explained a largely 11.3% and 15.4% of the phenotypic variance ( Table 3 ), which should be considered as major QTLs. About thirty QTLs of ONPO have been previously reported by linkage or association mapping in oilseed rape (Khan et al., 2019; Qadir et al., 2022; Ahmad et al., 2023). The results of comparative QTL analysis ( Supplementary Table S3 ) showed that four QTLs (qONPO.A01-1, qONPO.A10-1, qONPO.C02-1, qONPO.C05-2) were overlapped with the reported ONPO QTLs, while the remaining six should be novel loci. For example, the genomic region qONPO.A01-1 detected in the present study was overlapped with two association loci identified in our previous study (Qadir et al., 2022).

The previous studies have showed that four phytohormones (AUX, BR, CK, and GA) play a central role in regulating ONPO (Qadir et al., 2021). Therefore, the concentrations of these four phytohormones in ovaries at ovule initiation stage were measured and compared between high and low ONPO pools ( Table 4 ). The two pools respectively comprised ten high and nine low ONPO lines that were chosen from top or bottom 10% lines common in the investigated Spring and semi-Winter rapeseed environments, with an average difference of 12.6 and 9.4, respectively.

For auxin, the concentrations of IAM and ILA in H pool (3.52 and 4.90 pg/mg) were significantly lower than those in L pool (5.19 and 5.35 pg/mg). For brassinosteroid, the levels of 28-norBL and TY in H pool (0.08 and 0.20 pg/mg) were significantly higher compared to L pool (0.06 and 0.18 pg/mg). For cytokinin, the contents of DHZ7G and tZR in H pool (0.06 and 0.76 pg/mg) were respectively lower and higher than that in L pool (0.08 and 0.47 pg/mg). For gibberellin, the levels of GA1, GA9 and GA24 in H pool (0.43, 0.57 and 1.59 pg/mg) were significantly higher than those in L pool (0.33, 0.34 and 0.89 pg/mg). In conclusion, the concentrations of nine sub-types of AUX, BR, CK, and GA showed significant difference between H and L pools, which should be responsible for the ONPO difference in these lines.

In the reference genome of DarmorV4.1, a total of 289 genes were found to be the homologues of known ovule number genes, by using BLAST analysis ( Supplementary Table S8 ). Of these, 15 were located within the genomic regions of identified ONPO QTLs except for qONPO.A06-1 ( Table 5 ), which should be the underling candidate genes. In addition, the number of DEGs in each QTL region varied from 11 to 73, with the sum of 327 ( TSupplementary Table S5 ). More importantly, three DEGs in the QTL region (BnaA01g23050D, BnaA06g30720D and BnaC05g21200D) were also homologous to the known genes regulating ONPO, which should be considered as the most probable candidate genes prioritized for further functional validation. For example, BnaA01g23050D was located in the genomic region of qONPO.A01-1, its expression abundance in H pool was 3.76-fold that of L pool.

Previously, about 30 QTLs of ONPO ( Supplementary Table S3 ) have been reported by a few linkage or association mapping studies using SNP array in oilseed rape (Ahmad et al., 2023; Khan et al., 2019; Qadir et al., 2022). These ONPO QTLs were distributed in 14 of all the 19 linkage groups, which accounted for 1.22-17.38% of the phenotypic variance with additive effect ranging from -1.03 to 2.70. In this study, a total of 10 QTLs were identified in a DH population of 201 lines through linkage mapping across four environments, which were distributed on eight chromosomes, accounting for 7.0-15.9% of phenotypic variance ( Table 3 ). Of these, two QTLs in A07 and C02 linkage groups (qONPO.A07-1 and qONPO.C02-1) were repeatedly identified in different environments and showed a relatively large effect, which should be considered as major QTLs and important target for practical application and further gene cloning. The relatively low reproducibility of detected QTL was highly accordant with the relatively low heritability of ONPO due to larger environmental difference in comparison with previous studies (Ahmad et al., 2023; Khan et al., 2019). Furthermore, the comparative QTL analysis ( Supplementary Table S3 ) indicated that four QTLs in the present study were overlapped with those reported previously (Khan et al., 2019; Yuan and Kessler, 2019; Ahmad et al., 2023), which should be important targets for maker-assisted selection. More importantly, the remaining six QTLs should represent novel loci for ONPO. Except for qONPO.A07-1, the additive-effect direction of other five loci were all positive, indicating that favorable alleles were mostly from high-ONPO parent 3S1235. In addition, the high-ONPO genotype of peak SNP markers of these loci can be used for the molecular improvement of over number.

The ONPO is affected by both genetic environmental factors, including growth conditions, flower position (Qadir et al., 2022) and size (Wetzstein et al., 2013), and nutrient availability (Chalcoff and Aizen, 2016). Our recent review article summarized the reported ovule number genes and its regulatory pathways, which demonstrated that four types of phytohormones (AUX, CKs, GAs and BRs) played a central role (Qadir et al., 2021). It involves many biological process including phytohormones synthesis, degradation, transport, perception, signal transduction and downstream response etc., which regulate ovule development and finally affect its number. In our recent study (Qadir et al., 2022), the abundances of ABA, BA, GA4, IAA and JA showed significant difference between the two pools of high and low ONPO lines from an association population. To further reveal the roles of different sub-types of phytohormones in regulating ovule number, 48 sub-types of the above-mentioned four types of phytohormones were measured in the current study, of which nine had significant differences between two pools of high and low ONPO lines ( Table 4 ). In addition, a total of 159 DEGs between these two pools were classified into the functional category of phytohormone metabolism ( Supplementary Table S5 ), which involved eight types of phytohormones (ABA, AUX, BR, CK, ETH, GA, JA and SA). More importantly, of the total of 54 DEGs homologous to known genes regulating ovule number, at least 22 (40.7%) were relevant to phytohormones, of which 14 were classified into hormone metabolism and 8 belonged to transcription factors that contained AUX and ETH response element. These research results highly indicated that phytohormones should play a major role in regulating ONPO.

In the reference genome Darmor-bzh, 289 annotated genes ( Supplementary Table S8 ) were the homologues of 68 reported genes regulating ONPO (Qadir et al., 2021). Of these 289 homologues, 54 (18.69%) were differentially expressed between the two pools of high and low ONPO lines ( Supplementary Table S5 ). For an example, BnaA01g23050D is an orthologue of HAP13, and its down-regulation decreases ONPO (Wang et al., 2016). Among a total of 111479 annotated genes in the Darmor-bzh reference genome, only 7689 (6.86%) were differentially expressed between the two pools mentioned above ( Figure 4C ). The former proportion was nearly three folds of the latter, which suggested that these DEGs were more preferentially to ovule number genes. In addition, the reported ONPO genes were classified into eight functional categories (Qadir et al., 2022), including RNA (23), hormone metabolism (22), development (12), signalling (6), protein (4), miscellaneous (3), micro RNA (1) and tetrapyrrole synthesis (1). It should be noted that nearly half of these DEGs belonged to the abovementioned categories ( Supplementary Table S5 ), which suggested that they are highly involved in known pathways regulating ovule number. These results strongly supported that the DEGs identified in current study were tightly associated with ovule number.

Genetic variations are valuable resources for understanding the biological basis of important traits and developing novel varieties with favorable traits. In the current study, a total of 10 ONPO QTLs have been identified using a DH population that was genotyped by 50K SNP array and phenotyped in four environments. Through BLAST analysis with the reported genes regulating ONPO (Qadir et al., 2021), a total of 13 homologues were found in these QTL regions ( Table 5 ), which should be considered as the direct candidate genes for ONPO QTL after the confirmation of genetic sequence differences between two parents. It should be noted that eight of the 13 homologues were relevant to phytohormones, of which six were key enzymes in the metabolism of phytohormones and two were auxin or cytokinin response factors. This result was highly accordant with the observed phytohormone content difference between the two pools of high and low ONPP lines ( Table 4 ), which highlight the key roles of phytohormones in regulating ONPO in oilseed rape. In addition, a lot of studies have proved that the combination of QTL mapping and transcriptomic analysis is an efficient approach for screening candidate genes (Hussain et al., 2021). In the current research, a total of 327 DEGs were also located in the physical regions of detected ONPO QTLs, which should also be potential candidate genes for ONPO QTL ( Table 5 ; Supplementary Table S5 ). More importantly, a few of genes were common between DEGs and ovule number homologous genes, which should be considered as the most important candidate genes. The subsequent research should focus on the fine-mapping of major QTL followed by functional verification of the candidate genes.