Association Study of Agronomic and Oil Traits in Rapeseed (Brassica napus L.) Mutant Lines Using Genotyping-by-Sequencing

1. Introduction

Rapeseed (Brassica napus L.), an interspecific amphidiploid hybrid (2n = 38) [1], belonging to the genus Brassica in the family Brassicaceae, emerged from natural hybridization between Brassica rapa and Brassica oleracea [1,2]. This versatile crop is distinguished by its capacity to yield substantial amounts of oil, which is deemed a vital and salubrious vegetable oil worldwide [2]. Its multifarious applications span a wide range of industries, including but not limited to, food, animal feed, energy, and chemical industries [1,3]. As such, breeders have endeavored to breed new rapeseed cultivars with modified and advantageous agronomic characteristics, such as heightened yields and disease resistance, and oil traits, such as altered fatty acid compositions and crude fat content [2,4]. However, the genetic resources of rapeseed with natural properties that are germane to its agronomic and oil-related traits are meager in Korea [4]. Given the scarcity of genetic diversity in this crop, mutagenesis offers an efficacious approach to generate novel genetic resources, enabling the attainment of desirable traits, such as elevated yields and modified fatty acid compositions in rapeseed [1,2,4].

The weight of rapeseed seeds is widely considered to be a critical determinant of overall yield, and is thus a key trait for plant evolution and crop enhancement [4,5]. Seeds that are small in size have a greater potential for dispersion, while larger seeds tend to exhibit higher adaptability to both biotic and abiotic stress factors [2,6]. Additionally, the seedlings that emerge from larger seeds often possess superior competitive survival rates compared to those arising from smaller seeds, thus making seed weight a critical trait for plant fitness and survival [2,5,6]. Many crops have larger seed sizes than their wild relatives [5,6,7]. Hence, the identification of genes that are linked with thousand seed weight is of significant importance.

The quality of rapeseed oil is a crucial factor for the economic feasibility of its industrial application [3,4]. The nutritional and healthy properties of rapeseed oil are predominantly determined by the fatty acid compositions, which are synthesized through a series of biochemical reactions involving acetyl-CoA and NADPH [4,8]. Rapeseed oil, known for its high erucic acid content, was initially used in limited quantities due to its damaging effects on animal cardiac muscles [3,9,10]. However, modern breeding strategies for edible oil production have focused on developing canola-type rapeseed with double low seeds, low in erucic acid and glucosinolate content [1,10]. The resulting fatty acid composition of canola-type rapeseed oil, comprising 7% palmitic acid, 2% stearic acid, 61% oleic acid, 11% linoleic acid, and 21% linolenic acid, is considered nutritionally optimal, with high-oleic acid oil ( > 70%) gaining popularity as a source of healthy and stable cooking oil [1,10,11]. On the other hand, rapeseed oil with high erucic acid content has been found to be suitable for various industrial applications, such as polyethylene films, biodegradable plastic, biodiesel, printing, and the steel metal industries [10,12]. The crushing industry derives its primary value from the oil content of oilseed rapeseeds, despite the value of the protein for animal feed [3,10,11]. This is particularly crucial for the development of the biodiesel market and cooking oil production, where the optimization of production costs is of utmost importance [1,11,12]. Due to the tremendous economic importance of rapeseed oil and its expanding use as a renewable feedstock, the increased seed oil content is critical. Therefore, optimizing the fatty acid composition of rapeseed oil is a crucial objective of various breeding programs.

Mutation breeding is a method that involves the use of physical and/or chemical mutagens to induce genetic alterations in plants, resulting in the development of useful characteristics for commercialization purposes [4,13]. This process requires selecting the most effective mutation for modifying agronomic characteristics and oil traits, which can be achieved through various mutagenesis techniques [7,14]. Gamma-ray irradiation is one of the most commonly used techniques in plant mutation breeding due to its ability to induce vast amounts of genetic alterations [13,15]. It has been successful in inducing translocations in amphiploidy with desirable genes, resulting in seed yields, semi-dwarf varieties, and broader disease resistance [6,7,13]. Additionally, mutagenesis by gamma irradiation has been shown to stimulate genetic recombination, which broadens the mutation spectrum and enhances the efficacy of the technique [4,7,15].

Nowadays, advances in next-generation sequencing (NGS) technology have enabled effective sequencing of plant genomes and the direct detection of SNPs, which promotes the development of cultivars of interest in plant breeding [4,16]. GBS reduces the genomic complexity by digesting the genome into small fragments (with restriction enzymes) that are read on a short-read sequencing platform [16,17]. As the availability of entire genome sequences and SNP arrays has increased, association mapping has become a robust approach for determining genetic characteristics, which can significantly increase the precision of quantitative trait locus (QTL)-based position estimates [17,18]. Moreover, association mapping has been used to overcome the limitations of QTL mapping, particularly because of the large number of SNP markers identified by NGS [16,18]. The recent completion of the B. napus genome enabled direct com parisons between documented complex traits associated with mapping studies [4,7,11]. The identification of SNPs through GBS not only enables the analysis of genetic diversity but also facilitates the integration of GWAS into a unified research project [16,18,19]. GWAS are a powerful tool for identifying QTLs and genes in crops, such as rapeseed, soybean, and rice [6,17,19,20].

We have developed mutant rapeseed lines derived from gamma-radiation mutation, and these lines have a variety of agronomic characteristics (plant height, ear length, thousand seed weight and seed yield) and oil trait (fatty acid compositions and crude fat content). The objectives of this study were the identification of SNPs in 95 rapeseed mutant lines, and to identify candidate genes related to agronomic and oil traits using a GWAS.

2. Materials and Methods

2.1. Plant Material

The seed of the ‘Tamra’ cultivar was obtained from the Bioenergy Crop Research Center (Rural Development Administration, Korea). Mutant rapeseed lines were generated by treating the seeds with 700 Gy of gamma (⁶⁰Co) radiation at 2008 [4]. The procedure used to develop mutant rapeseed lines is shown in Figure 1. The treated seeds were sown to obtain the M₁ generation, and seeds from one silique (developed from the main stem of each M₁ plant) were harvested. M₂ seeds from 200 individual plants were grown with a single replicate. In the M₂ generation, all individuals were investigated for morphological and agronomic mutations relative to the original cultivar. One hundred and twenty rapeseed mutants, selected based on their agronomic characteristics were obtained from the M₃ and M₅ generations. We analysed the uniformity of the fatty acid compositions by GC-MS (Gas chromatography mass spectrometry) and crude fat content for two generations (M₆ to M₇) to select stable lines. Finally, 95 rapeseed mutants lines that varied in agronomic characteristics, fatty acid compositions and crude fat content exhibited stable inheritance of the mutated characteristics from M₈ generations were selected. The selfing procedure was continued until the M₈ generation. Four agronomic traits including plant height, ear length, thousand seed weight and seed yield were investigated according to the International union for the protection of new varieties of plants (UPOV) test guidelines for rapeseed and the standard of research and investigation for agronomic traits. Cultivars with radiation-generated mutant genotypes were grown by the Radiation Breeding Research Team at the Advanced Radiation Technology Institute of the Korea Atomic Energy Research Institute, Korea.

2.2. Determination of Fatty Acid Compositions and Crude Fat Contents

The seed oil content was analyzed using the AOAC method as previously described. Using the Soxhlet extraction procedure, 5 g crushed seeds (80 mashed) were packed into a thimble and the oils were extracted with diethyl ether for 6 h. Fatty acid compositions were measured previously described method. The rapeseed oil was extracted from rapeseed powder in 1 mL of chloroform-hexane-methanol (8:5:2, v/v/v) for 12 h. From this, 200 μL of extracted oil was added to 75 μL of methylation reagent (0.25 M methanolic sodium methoxide: petroleum ether: ethyl ether, 1:5:2, v/v/v) for derivatization. Hexane was added to bring the total volume up to 1 mL. The fatty acid composition of the rapeseed seed oil was analyzed using a GC-MS (Plus-2010, Shimadzu, Kyoto, Japan) instrument equipped with an HP-88 capillary column (J&W Scientific, 60 m × 0.25 mm × 0.25 μm) under the following conditions: ionization voltage, 70 eV; mass scan range, 50–450 mass units; injector temperature, 230 °C; detector temperature, 230 °C; injection volume, 1 μL; split ratio, 1:30; carrier gas, helium; and flow rate, 1.7 mL/min. The column temperature program specified an isothermal temperature of 40 °C for 5 min increasing to 180 °C at a rate of 5 °C/min then a subsequent increase to 230 °C at a rate of 1 °C/min. It was identified that the substances present in the extracts according to their retention time (RT) and using the mass spectra database (NIST 62 Library).

2.3. DNA Extraction

Young leaves were sampled from the original cultivar ‘Tamra’ and 95 rapeseed mutant lines. Genomic DNA was isolated using a DNeasy 250 Plant Mini Kit (Qiagen, Hilden, Germany) according to the manufacturer’s instructions. The extracted DNAs were stored at 4 °C until use. For polymerase chain reaction (PCR) analysis and, DNA concentrations were determined using a NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific., Waltham, MA, USA) were then adjusted to 30 ng/μL.

2.4. Library Construction and Genotyping-by-Sequencing (GBS)

GBS libraries were constructed using the restriction enzyme ApeKI (5′-GCWGC-3′; New England Biolabs, MA, USA) following a protocol modified from that of a previous study [21]. Oligonucleotides constituting the top and bottom strands of each barcode adapter and a common adapter were separately diluted with TE buffer (50 μM each) and annealed using a thermocycler. DNA samples (100 ng/μL) were added to individual adapter-containing wells. Samples (DNA with adapters) were digested with ApeKI overnight at 75 °C. Sets of digested DNA samples, each with a different barcode adapter, were combined (5 μL each) and purified using the QIA quick PCR Purification Kit; (Qiagen, CA, USA), according to the manufacturer’s instructions. Restriction fragments from each library were then amplified in 50 μL volumes containing 2 μL of pooled DNA fragments, HerculaseII Fusion DNA Polymerase (Agilent, CA, USA), and 25 pmol each of the following primers: (A) 5′-AAT GAT ACG GCG ACC ACC GAG ATC TAC ACT CTT TCC CTA CAC GAC GCT CTT CCG ATC T-3′ and (B) 5′-CAA GCA GAA GAC GGC ATA CGA GAT CGG TCT CGG CAT TCC TGC TGA ACC GCT CTT CCG ATC T-3′. These amplified sample pools constituted a sequencing ‘library’. The library was sequenced on the Illumina Hiseq X ten platform by SEEDERS Co. (Daejeon, Korea).

2.5. Sequence Preprocessing and Alignment to Reference Genome Sequence

Demultiplexing was performed using the barcode sequence, and adapter sequence removal and sequence quality trimming were performed. Adapter trimming was performed using Cutadapt (ver. 1.8.3) [22], and sequence quality trimming was performed using DynamicTrim and LengthSort of the SolexaQA program (ver.1.13) [23]. DynamicTrim cuts low-quality bases at both ends of short reads according to the Phred score, and refines it with high-quality cleaned reads. LengthSort removes excess base cuts made in DynamicTrim; Phred score of Dynamic-Trim ≥ 20, and LengthSort using short read lengths ≥ 25 bp. BWA (0.7.17-r1188) [24] generated cleaned reads passing the preprocessing process, and performed mapping to the reference genome sequence. Mapping was a preliminary step to detect raw SNPs (In/Del) sequenced samples. BAM format file was created, and default values were used, except for the following options: a seed length (−l) of 30, maximum differences in the seed (−k) of 1, number of threads (−t) of 16, mismatch penalty (−M) of 6, gap opening penalty (−O) of 15, and gap extension penalty (−E) of 8. The experiment was conducted in repetition.

2.6. Raw SNP Detection and Consensus Sequence Extraction

Clean reads were mapped to the standard genome sequence, and the generated BAM format files were used to detect raw SNPs using SAMtools (ver. 0.1.16) [25], and extract consensus sequences. SNP validation was performed using SEEDERS in-house script [26] before SNP detection; raw SNP detection was performed, and default values were used, except for the following options: a minimum mapping quality for SNPs (−Q) of 30, minimum mapping quality for gaps (−q) of 15, minimum read depth (−d) of 3, minimum indel score for nearby SNP filtering (−G) of 30, SNPs within INT bp around a gap to be filtered (−w) of 15, window size for filtering dense SNPs (−W) of 15, and maximum read depth (−D) of 675.

2.7. Generate SNP Matrix

To conduct the analysis of SNPs between the analyzed objects, an integrated SNP matrix was generated between samples. A list of unions was created using the raw SNP positions obtained by comparing each sample with a standard dielectric, and a non-SNP locus was filled in from the consensus sequence of the sample. Then, the final SNP matrix was generated by filtering the mis-called SNP positions through SNP comparison among samples. SNPs were divided into homozygous (SNP read depth ≥ 90%), heterozygous (40% ≤ SNP read depth ≤ 60%), and other (homozygous/heterozygous; could not be distinguished by type) groups based on their position. The named SNP positions were classified into “intergenic or genic regions” based on the position information of the standard genome sequence, and the genic region was further classified into “CDS or intron regions”. The common SNP in the original cultivar ‘Tamra’ was selected first in the integrated SNP matrix position between mutant lines, and the polymorphic SNP was selected by comparing the common SNP of the original cultivar with the base sequence of each mutant. In order to perform gene ontology (GO) analysis and flexible relationship analysis, SNP loci of each mutant line were integrated to secure the SNP locus of the union.

2.8. Gene Ontology (GO) Analysis of Genes with Polymorphic SNPs

Gene ontology alignment was performed using candidate sequences containing polymorphic SNPs and sequences provided by GO database in house scripts [27]. Thresholds were classified into three functional categories: BP (Biological Process), CC (Cellular Component), and MF (Molecular Function), with a significance level of 0.01 (E-value ≤ 1.0 × 10⁻¹⁰, best hits).

2.9. GWAS with Agronomic Characteristics, Fatty Acid and Crude Fat

For the GWAS, a total 73,226 filtered SNPs with 5% < minor allele frequency, missing data < 30% were filtered from raw data and were used to perform GWAS analysis using generalized linear model (GLM) in Trait Analysis by Association, Evolution and Linkage (TASSEL), TASSEL 5 [28]. All parameters for GWAS were set as default values. The –log₁₀(p) threshold of the quantile-quantile plot and Manhattan plot was calculated using the Bonferroni method (p = α/n). With 73,226 SNPs employed in this study, at α = 1, the Bonferroni-corrected thresholds for the p values were 1.36 × 10⁻⁵ (α = 1), with corresponding –log₁₀(p) values of 4.864 for the suggestive threshold [29].

3. Results

3.1. Agronomic Characteristics

Ninety-five rapeseed mutant lines and their original cultivar ‘Tamra’ were selected for evaluation in terms of their agronomic characteristics (Table 1). The rapeseed mutant population showed wide variability in the traits measured over the course of 2 years. The plant height of the original cultivar was 163.5 cm. The plant height of rapeseed mutant lines ranged from 134.0 cm (Tr2-2) to 175.0 cm (Tr25-14). The ear length of the original cultivar was 55.5 cm. The ear length of all mutant lines ranged from 35.0 cm (Tr8-3-1) to 77.5 cm (Tr38-7). Significant differences in the thousand seed weight and the seed yield were observed between all mutant lines. With the original cultivar, the thousand seed weight was 3.8 g. The Tr6-11-1 line had the highest (5.4 g), whereas the Tr138-L line had the lowest (2.8 g) thousand seed weight. The seed yield of original cultivar was 309 kg/10a. The highest seed yield (398 kg/10a) was observed with the Tr38-4 line, and the lowest seed yield (144 kg/10a) was found with the Tr14-3 line.

3.2. Fatty Acid Compositions and Crude Fat Contents

The fatty acid compositions and crude fat content of the original cultivar ‘Tamra’ and 95 rapeseed mutant lines genotypes are shown in Table 2. The composition of palmitic (C16:0), stearic (C18:0), oleic (C18:1), linoleic (C18:2), linolenic (C18:3), and eicosenoic acid (C20:1) of the original cultivar was 5.5%, 1.9%, 65.5%, 20.5%, 5.6%, and 1.0%, respectively. Palmitoleic acid (C16:1) and erucic acids (C22:1) were not detected in the original cultivar. On the other hand, 0.1% to 0.76% of palmitoleic acid was detected in 50 mutant lines. The stearic acid content of the mutant lines ranged from 0.2% to 4.4%, with an average of 1.7%. Significant differences in oleic acid composition were observed among the mutant lines. The oleic acid composition of the mutant lines ranged from 32.4% to 81.7%, and the oleic acid content was the highest in the Tr4 mutant line. Likewise, nine mutant lines (Tr2-late, Tr3, Tr5, Tr5-1, Tr7, Tr25-10, Tr8-3-1, Tr6-4, and Tr8-1) also showed higher oleic acid composition ( > 70%) than the other genotypes. The highest linoleic acid composition was recorded at 33.5% in the Tr138-15 mutant line, and the lowest value of 7.9% was found in the Tr1-1 line. The highest linolenic acid content was found in the Tr5-3 line (13.7%), whereas the lowest content was detected in the Tr65-7 line (4.0%). Eicosenoic acid was not present in seven mutant lines (Tr25-12, Tr38-7, Tr38-16, Tr38-17, Tr65-7, Tr65-9, and Tr65-10). The eicosenoic acid content of the remaining mutant lines ranged from 0.4% to 13.1%, with the highest composition recorded in Tr14-1. Erucic acid was detected in 43 mutant lines with a composition range of 0.2% to 29.0%. Two mutant lines (Tr14-1 and Tr14-5) showed a relatively high composition of erucic acid ( > 20%). The crude fat content of the original cultivar was 40.16 mg/100 g. The crude fat content for all mutant lines ranged from 21.31 to 40.98 mg/100 g with an average of 37.22 mg/100 g. Crude fat content similar to that of the original ‘Tamra’ cultivar was observed in 9 mutant lines (Tr5-1, Tr14-8, Tr14-9, Tr14-20, Tr25-6, Tr25-10, Tr65-10, Tr138-13, and Tr8-5-1). The Tr14-5, Tr18-2, Tr25-1, and Tr25-2 lines contained lower levels of crude fat than the original cultivar.

3.3. GBS Analysis of Rapeseed Mutant Lines

A GBS library consisting of 96 rapeseed genotypes (95 novel mutant lines and the original cultivar) was sequenced using the Illumina HiSeq X ten platform (Illumina, Madison, WI, USA). A summary of GBS sequencing results is presented in Table 3 and Table S1. In total, 715 million reads comprising 108,004,241,578 nucleotides (108 Gb) were generated, with 7.45 million reads per genotype. After removing low-quality sequences, 655,243,166 clean reads remained, with 6.8 million reads per genotype on average. The total length of the clean reads ranged between 73,633,121 base pairs (bp) and 5,305,962,156 bp, with an average read length of 777,311,365 bp (Table S1). The total number of mapped reads was 651,288,040 for all lines, with an average of 6,784,250 reads per sample. The mapped-read rates (%) ranged from 98.98% to 99.46%. On average, 99.39% of the filtered reads were mapped to the reference genome sequence (Table S1). The total length of the mapped region was 3,163,629,539 bp, with an average of 32,954,474 bp per sample, which covered approximately 3.57% of the reference genome sequence.

3.4. Identification of SNPs

The SNPs for each genotype were first selected from the union of SNPs in the matrix position between the samples and the reference genome sequence. A summary of SNPs detection of the original cultivar ‘Tamra’ and 95 rapeseed mutant lines are shown Table S2. In total, 9,381,607 SNP loci were identified for all genotypes, of which 5,392,272 SNP loci were homozygous for the major allele, 1,500,793 SNPs were heterozygous, and SNP 2,448,542 SNPs were etc. (homozygous/heterozygous; could not be distinguished by type) (Table S2). The largest number of SNPs was recorded for the original cultivar (183,638 SNPs), and the lowest number was recorded for the Tr65-2 mutant line (17,097 SNPs). The number of homozygous SNPs ranged from 11,879 to 101,493, with an average of 56,170 per sample. The number of heterozygous SNPs ranged from 2,176 to 39,895, with an average of 15,633 per sample. In other words, more homozygous SNPs were found than heterozygous SNPs. The SNP distribution in the rapeseed chromosomes ranged from 4,870 (chromosome A04) to 14,497 (chromosome C03), with a mean of 9,378 per chromosome (Figure 2). Specifically, the number of SNPs was highest in the order of chromosomes A09, C01, C02, C03, and A09.

To identify gamma irradiation-induced mutations, the polymorphic SNPs were identified by comparing common SNPs (the first selected 9,381,607 SNPs) in the original cultivar with the corresponding sequences of the mutant lines (Table S3). In total, 182,257 union SNPs were observed among all mutant lines, where the largest number of polymorphic SNPs was recorded for Tr6-11-1 (25,568 SNPs), and the lowest number was recorded for the Tr65-2 mutant line (2,281 SNPs). The number of homozygous polymorphic SNPs ranged from 1,620 to 19,590, with an average of 8,225 per mutant line. The number of heterozygous polymorphic SNPs ranged from 661 to 14,277, with an average of 4,051 per mutant line.

Functional annotation of the SNPs in the rapeseed mutant lines was performed using the reference genome sequence (Table S4). In total, 181,225 union SNPs were classified for genome annotation. Of the 181,225 union SNPs, 99,098 SNPs (54.7%) were located in genic regions, and 82,127 SNPs (45.3%) were detected in intergenic regions, encompassing 30,121 genes. Within the genic regions, some SNPs were found in untranslated regions (UTRs, 3,166 SNPs), and the remaining SNPs were located in coding sequences (67,561 SNPs) and introns (29,067 SNPs).

3.5. Gene Ontology (GO) Analysis of Genes with Polymorphic SNPs

To functionally classify the genes mutated by gamma-ray irradiations in the mutant lines, GO enrichment analysis was conducted with 30,121 genes carrying polymorphic SNPs (p < 0.05) (Table S5). The genes were classified into three main categories: biological process (BP), cellular component (CC), and molecular function (MF) genes (Figure 3). Genes containing BP SNPs included those involved in cellular processes (13,144 genes), metabolic processes (11,154 genes), organic substance metabolic processes (10,584 genes), and cellular metabolic processes (10,038 genes). CC SNPs were found in genes associated with intracellular anatomical structures (15,349 genes), organelle entities (13,635 genes), and intracellular organelles (13,624 genes). MF SNPs were identified in genes related to GO terms such as binding (9,513 genes), organic cyclic compound binding (5,795 genes), heterocyclic compound binding (5,776 genes), ion binding (4,955 genes), and protein binding (3,451 genes).

3.6. GWAS for Agronomic Characteristics

In this study, association analysis was performed with a generalized linear model for the four agronomic traits. According to the results of the Manhattan and QQ plots presented in Figure 4, 76 SNPs across 14 chromosomes for TSW were detected at the significance threshold level (-log₁₀(p) = 4.864) (Table 4 and Table S6). However, no suggestive or significant SNPs were detected for the PH, EL, and SY. Among the 73,226 union SNPs dataset associations identified, 76 SNPs were significantly associated with agronomic characteristics (Table S7), when applying the generalized linear model (GLM). Of the 76 selected union SNPs, 55 were located in genic regions, and 21 were detected in intergenic regions. A total of 14 genes (S-formylglutathione hydrolase-like (SFGH), zinc finger BED domain-containing protein RICESLEEPER 1-like (ZBED1), midasin-like (MDN1), villin-3-like isoform X2 (VLN3), death-associated inhibitor of apoptosis 1 (DIAP1), ubiquitin-like-specific protease 1C (ULP1C), chitinase 1 (CHIT1), amino acid transporter AVT1A (AVT1A), endoribonuclease Dicer homolog 1-like (DLC1), AP-1 complex subunit mu-2 (AP1M2), RNA-binding protein 5 (RBP5), ATPase 6, plasma membrane-type (AHA6), serine acetyltransferase 5 (SAT5), and U2 snRNP-associated SURP motif-containing protein isoform X2 (U2SURP)) were annotated from thousand seed weight.

3.7. GWAS for Fatty Acid and Crude Fat

In this study, association analysis was performed with a generalized linear model for the 8 fatty acid and crude fat. According to the results of the Manhattan and QQ plots presented in Figure 5, 2 SNPs across 1 chromosome for C16:1, 1 SNP across 1 chromosome for C18:0, 1 SNP across 1 chromosome for the C18:1, 6 SNP across 6 chromosomes for C20:1, 26 SNPs across 10 chromosomes for C22:1, 37 SNPs across 10 chromosomes for crude fat were detected at the significance threshold level (-log₁₀(p) = 4.864) (Table 5 and Table S8). However, no suggestive or significant SNPs were detected for the C16:0, C18:2, and C18:3. Among the 73,226 union SNPs dataset associations identified, 73 SNPs were significantly associated with oil traits (Table S9), when applying the generalized linear model (GLM). Of the 73 selected union SNPs, 37 were located in genic regions, and 36 were detected in intergenic regions. One gene (Encodes a putative E3 ubiquitin ligase (SUD1)) was annotated from C16:1. For C20:1, two SNPs were located in the genic region (CDS), and were associated with probable Xaa-Pro aminopeptidase P isoform X2 (XPNPEP1), and hypothetical protein. A total of 26 SNPs (CDS; 6 SNPs, introns; 3 SNPs, intergenic; 16 SNPs, and UTR; 1 SNP) loci were significantly associated with C22:1, including six genes (EXO70 gene family (EXO70), paramyosin-like (Prm), mediator of RNA polymerase II transcription subunit 34 (MED34), cyclin-C1-1 (CCNC), ribonuclease T2 family (RNASET2), and protein ENHANCED DISEASE RESISTANCE 2-like isoform X2 (EDR2)) that were annotated from C22:1. A total of 10 genes (Part of a nanodomain complex that tethers PI4Kα1 to the plasma membrane (PI4Kα1), 1-aminocyclopropane-1-carboxylate synthase 6 (ACS6), mitoferrin (Mfrn), Transmembrane amino acid transporter family protein (TAAT), chromo domain-containing protein LHP1-like (LHP1), non-lysosomal glucosylceramidase (GBA2), RING-H2 finger protein ATL47-like (ATL47), embryogenesis-associated protein EMB8 (EMB8), CEP DOWNSTREAM 1 (CEPD1), and protein CHUP1, chloroplastic-like (CHUP1)) were annotated from crude fat.

4. Discussion

High seed yield has been most important factor in rapeseed breeding. Normally, seed yield revealed to the environmental conditions, genotypes and those interactions [5]. In a previous study, genotypes have a great influence on the rapeseed seed yield [5,6]. Also, the relatively limited genetic resources of rapeseed breeding materials restrict the potential for further improvement in seed yield [6,7,30]. The rapeseed seed yield depends on various agronomic traits, such as the ear length, branch numbers, silique numbers, number of seeds/siliques, thousand seed weight, and disease resistance [2,5,6,9,30]. The original cultivar ‘Tamra’ has low seed yield (300–350 kg/10a-1), resulting in low industrial value [31]. In this study, the original cultivar ‘Tamra’ and 95 rapeseed mutant lines were analyzed during two years for the assessment of yield. Fifteen rapeseed mutant lines were approximately 26% higher (over 390 kg/10a) than that of the original cultivar. These results indicate that it is useful as a material for developing new rapeseed cultivar with high seed yield trait.

Industrial demands for rapeseed oils includes the ability to change the quality of fatty acids, and the demands are recognized by various breeders [3,12]. Their values and major applications depend on their fatty acid compositions [3,4]. Rapeseed oil is rich in unsaturated fatty acids, such as oleic acid, linoleic acid, and linolenic acid, which are generally recognized as useful oil compounds [3,8,11,12]. Gamma rays have been widely applied to mutation breeding programs because of their high radiolucency to plant organs and can cause mutations affecting agronomic and chemical traits [13,14,32]. Furthermore, gamma rays are particularly effective in changing the oil traits of various oilseed crops [20,32]. Radiation-induced mutagenesis techniques have been successfully developed with oil crops to generate useful mutants with economic traits, oil contents, and fatty acid compositions [11,13,33]. Rapeseed has been mutagenized using physical mutagens to improve many important characteristics, including its ability to produce high-quality industrial oil [4,12,34]. Oleic acid (C18:1) and linoleic acid (C18:2) are considered healthy and nutritious [3,8]. On the other hand, the three double-bonds of extracted linolenic acid (C18:3) are easily oxidized, which leads to a reduced frying thermal stability and storage time of the oil [3,11]. It is also important to develop cultivars with a high percentage of erucic acid oil, a valuable raw material used in the plastic film and nylon industries [3,10,11]. The original cultivar ‘Tamra’ oleic acid content was 65.5%, the erucic acid content was 0%. The oleic acid contents of the Tr3 (78.6%), Tr4 (81.73%), Tr5 (79.43%), and Tr7 (78.7%) mutant lines were higher ( ≥ 75%) than those of different genotypes. The erucic acid contents of the Tr14-1 (23.6%) and Tr14-5 (29%) mutant lines were higher ( ≥ 20%) than those of different genotypes. In this study, oleic acid and erucic acid showed the greatest changes. Tr14-1 had an erucic acid content of 23.6% when the oleic acid content was 32.4%, and Tr14-5 had an erucic acid content of 29% when the oleic acid content was 32.8%. Thus, the gamma ray irradiation induced mutation influence on fatty acid compositions in rapeseed [3,4,13,20] and these mutant lines are potentially useful as materials for developing new rapeseed cultivars with improved oil traits.

GO analysis of rapeseed genes containing polymorphic SNPs induced by gamma irradiation was conducted. GO grouping of the mutated genes indicated that many variations induced functional changes. With respect to the BP category, mutations occurred most frequently in genes involved in cellular processes (GO:0009987), metabolic processes (GO:0008152), organic substance metabolic processes (GO:0071704), cellular metabolic processes (GO:0044237), and primary metabolic processes (GO:0044238). Regarding the CC category, the associated intracellular entities included intracellular anatomical structure (GO:0005622), organelle (GO:0043226), intracellular organelle (GO:0043229), membrane-bound organelle (GO:0043227), and intracellular membrane-bound organelle (GO:0043231). In terms of the MF category, the major polymorphic SNPs were linked to GO terms such as binding (GO:0005488), organic cyclic compound binding (GO:0097159), heterocyclic compound binding (GO:1901363), ion binding (GO:0043167), and protein binding (GO:0005515). Guan et al. (2023) reported that GO enrichment analysis of genes with high expression levels during seed development showed in B. napus that the top three terms were gene expression, translational initiation, and cellular nitrogen compound metabolic process in the BP category; intracellular anatomical structure, organelle, and intracellular organelle in the CC category; and translation regulator activity, translation factor activity, RNA binding, and translation regulator activity, nucleic acid binding in the MF category [35]. These results suggest that the top three GO terms may play an important role in seed development in B. napus, and our results confirmed all of these terms (Table S5). Previously, GO analysis of variation in fatty acid-mutated 46 rapeseed genotypes derived from gamma rays included cellular processes (BP), primary metabolic processes (BP), and nitrogen compound metabolic processes (BP); intracellular entities (CC) organelles (CC) and intracellular organelles (CC); and nucleotide binding (MF), nucleoside phosphate binding (MF), and anion binding (MF) [4]. Therefore, it appears reasonable to expect that an appreciable similarity in mutations is observed in radiation mutagenesis among rapeseed mutant lines.

The GWAS detected significantly (-log₁₀(p) > 4.864) associated 14 SNPs with thousand seed weight. Seed weight is a complex trait that is regulated by multiple genes and influenced by genotypes [5,6]. Seed weight is an important agronomic trait as it affects seed yield, germination rate, and seedling vigor [5,6,7,10]. Therefore, understanding the molecular mechanisms underlying seed weight regulation is of great interest to plant breeders. Several genes have been identified that play a role in controlling seed size and weight in plants process. Our findings suggest that SFGH, ZBED1, MDN1, VLN3, DCL1, AP1M2, RBP5, AHA6, SAT5, AVT1A, U2SURP, DIAP1, ULP1C, and CHIT1 are associated with thousand seed weight. The SFGH is a gene that encodes a protein that belongs to the esterase/lipase/thioesterase family [36]. This gene is known to play a role in the detoxification of formaldehyde, which is no direct evidence linking SFGH to seed weight, there have been studies that suggest that SFGH may play a role in larger leaves, higher biomass, and increased grain yield in wheat [36,37]. The relationship between SFGH and rapeseed seed weight is not well-understood and further research is needed to fully elucidate the role of this gene. The ZBED1 gene is a transcription factor that plays an important role in regulating seed development and maturation that overexpression of ZBED1 in rice (Oryza sativa) resulted in increased seed weight and yield [38]. Studies have shown that ZBED1 regulates the expression of genes involved in lipid metabolism, hormone signalling, and stress responses, which are all important processes in seed development and maturation [38,39]. Therefore, it is possible that the ZBED1 gene could have an impact on seed weight. The MDN1 gene as “low quality protein” genes are believed to play a role in various cellular processes, such as protein quality control and regulation of gene expression [40]. MDN1 has been shown to play a role in regulating seed weight by affecting the accumulation of storage proteins in the seeds in Arabidopsis thaliana [41]. Similarly, in rice, a MDN1 gene has been shown to affect seed weight by regulating the expression of genes involved in the synthesis of storage proteins [40,41]. VLN3 are actin-binding proteins that play important roles in cytoskeletal organization, cell migration, and cell signalling in various organisms, including plants [42]. VLN3 gene have been found to be associated with seed development and germination [42,43]. Therefore, VLN3 may play a role in seed development, more research is needed to determine its specific functions and the extent of its involvement in seed weight regulation. The DCL1 is involved in the process of RNA interference (RNAi), which is a mechanism for gene regulation [44]. Several studies have suggested a relationship between DCL1 and seed weight that loss of function mutations in DCL1 result in smaller seed size and low weight, while overexpression of DCL1 leads to larger seeds in Arabidopsis thaliana and rice [45]. One proposed mechanism for the role of DCL1 in seed weight regulation is through its involvement in the processing of microRNAs (miRNAs) [44,45]. The AP1M2 is a transcription factor that plays a role in various biological processes, including cell differentiation, proliferation, and apoptosis and the AP1M2 gene complex has also been linked to seed development and seed weight regulation [46]. The AP1M2 that has been implicated in seed weight regulation is the JUN-like transcription factor (JUB1) [47]. JUB1 has been shown to promote seed size and weight by regulating the expression of genes involved in cell proliferation and expansion in Arabidopsis thaliana [46,47]. Additionally, the AP1M2 has been shown to regulate seed size by promoting cell division and elongation in soybean [46,47]. The RBP5 is a protein that binds to RNA molecules and plays a role in post-transcriptional regulation of gene expression [48]. Previous studies have suggested that RBP5 may directly affect seed weight by overexpression of RBP5 resulted in larger seeds, likely due to increased expression of genes involved in embryo development and carbon metabolism in rice and Arabidopsis [48,49]. These studies suggest that RBP5 may play a role in this process by regulating the expression of genes involved in seed development and lipid metabolism. The AHA6 gene is a protein that belongs to the P-type ATPase family, which is involved in ion transport across membranes [50]. AHA6 has been studied in relation to seed weight regulation that overexpression of AHA6 resulted in increased seed weight and size, likely due to increased nutrient uptake and transport to the developing seeds [51]. Additionally, AHA6 was shown to interact with a protein involved in seed development and to affect the expression of genes involved in lipid metabolism and seed coat development in Arabidopsis [50,51]. These findings suggest that AHA6 may play a role in seed weight regulation by facilitating nutrient uptake and transport to developing seeds and by influencing the expression of genes involved in seed development and metabolism [50,51]. Therefore, AHA6 affects seed weight and to determine its potential as a target for crop improvement. SAT5 activity is positively correlated with seed weight that overexpression of SAT5 resulted increased to seed weight, while reducing SAT5 activity led to a decreased seed weight in maize [52]. The AVT1A gene encodes an amino acid transporter, which is involved in the uptake and transport of amino acids [53]. As for the relationship between AVT1A and seed weight, some studies have shown that mutations in the AVT1A gene can affect seed weight in Arabidopsis thaliana and Brassica napus [53,54]. These studies found that mutations in the AVT1A gene resulted in smaller seeds with reduced amino acid content [53,54]. This suggests that AVT1A plays a role in regulating seed size and amino acid metabolism in plants. The U2SURP gene is a protein involved in pre-mRNA splicing in plants [52]. As for the relationship between U2SURP and seed weight, a study on Arabidopsis thaliana found that overexpression of U2SURP resulted in early flowering and reduced root growth, which in turn led to reduced seed weight [55]. This suggests that U2SURP may play a role in regulating seed weight indirectly by influencing other aspects of plant growth and development [52,55]. The DIAP1, ULP1C, CHIT1 gene are currently no known direct relationship with seed weight. This concluded that SFGH, ZBED1, MDN1, VLN3, DCL1, AP1M2, RBP5, AHA6, SAT5, AVT1A, and U2SURP may be associated with the thousand seed weight.

The GWAS detected significantly (-log₁₀(p) > 4.864) associated 8 SNPs with three fatty acid (palmitoleic acid, eicosenoic acid and erucic acid). The encodes a putative SUD1 showed significant associations with palmitoleic acid. However, the palmitoleic acid in relation to SUD1 is not clear. E3 ubiquitin ligases are enzymes that are involved in the regulation of protein degradation and play a key role in various cellular processes [56]. While, previous studies have suggested that palmitoleic acid may be involved in the regulation of ubiquitin-mediated protein degradation, more research is needed to understand the mechanisms underlying this process [56,57]. One gene XPNPEP1 was associated with the eicosenoic acid, these gene is no direct connection eicosenoic acid and fatty acid biosynthesis. A total of six genes (EXO70, Prm, MED34, CCNC, RNASET2 and EDR2) were annotated from erucic acid. The EXO70 encodes for proteins that are involved in the formation of the exocyst complex, which is responsible for the targeted delivery of membrane vesicles to specific locations on the plasma membrane [58]. The exocyst complex is essential for a variety of cellular processes, including cell polarity, cell division, and secretion [59]. The EXO70 encodes for a group of proteins that are involved in the formation of the exocyst complex, which is essential for the transport of secretory vesicles to the plasma membrane in eukaryotic cells [58,59]. While the role of the EXO70 family in exocytosis is well established, there is also evidence suggesting that some members of the family may be involved in lipid transport and secretion metabolism, the expression of EXO70 family (EXO70H4) was upregulated in response to glucose and that plants with mutations in the gene had reduced levels of fatty acids in Arabidopsis thaliana [58,60]. Another study found that the EXO70A1 protein in the filamentous fungus Aspergillus nidulans is involved in the regulation of lipid metabolism that mutations in the gene resulted in altered lipid composition, including a decrease in fatty acid content [58,59,60]. The specific mechanisms by which the EXO70 family members regulate rapeseed fatty acid biosynthesis are not yet fully understood, these studies suggest that they may play a role in this important metabolic process. MED34, Prm, CCNC and RNASET2 are no known direct relationship with erucic acid. EDR2 gene plays a role in regulating the plant immune system, helping it to defend against pathogen attacks [61]. EDR2 has been shown to be involved in regulating fatty acid biosynthesis in Arabidopsis thaliana, where mutations in the EDR2 gene led to increased accumulation of unsaturated fatty acids in plant tissues [62]. EDR2 acts as a negative regulator of fatty acid biosynthesis by repressing the expression of several genes involved in this process [61,62]. It is thought that EDR2 does this by interacting with other proteins involved in lipid metabolism and signaling pathways [61,62]. Further research is needed to fully understand its molecular mechanisms and potential applications in rapeseed. A total of 10 genes (PI4Kα1, ACS6, Mfrn, TAAT, LHP1, GBA2, ATL47, EMB8, CEP D1, and CHUP1) were annotated from crude fat content. PI4Kα1 is catalyzes the phosphorylation of phosphatidylinositol to produce phosphati-dylinositol 4-phosphate (PI4P) [63]. PI4Kα1 is an important signalling lipid that regulates various cellular processes, including vesicular trafficking, signal transduction, and cytoskeleton organization [63,64]. GBA2 is involved in the breakdown of a type of lipid molecule called glucosylceramide. This gene is not located in the lysosomes, which are specialized cellular compartments involved in the degradation of various macromolecules [65]. ATL47 is a type of E3 ubiquitin ligase, which has been found to play a role in regulating plant growth and development, as well as in response to regulating the activity of enzymes involved in fatty acid metabolism, such as acyl-CoA synthetases and acyl-CoA oxidases [66,67]. Moreover, some studies have reported that the overexpression of RING-H2 finger proteins, including ATL47, can lead to an increase in lipid accumulation in plant cells, suggesting a potential role in promoting fat biosynthesis [66,67]. The CEPD1 fat biosynthesis refers to the process of producing fatty acids and other types of lipids in living organisms. This process occurs downstream of the CEP (Carnitine/acylcarnitine transporter, EF-hand protein) pathway, which plays an important role in the transport of fatty acids into the mitochondria for energy production [68]. In the downstream process of fat biosynthesis, acetyl-CoA molecules are combined to form long-chain fatty acids through a series of enzymatic reactions [4,8,9,10]. ADP-glucose pyrophosphorylase is not directly involved in lipid biosynthesis. However, there is some evidence to suggest that there may be a connection between carbohydrate metabolism and lipid biosynthesis as excess glucose can be converted to acetyl-CoA, which is a precursor for fatty acid synthesis, through a process called de novo lipogenesis [9,11,12]. In addition, ADP-glucose pyrophosphorylase activity has been shown to be regulated by the availability of acetyl-CoA, which is an intermediate in fatty acid biosynthesis. It has been proposed that changes in acetyl-CoA levels may modulate the activity of ADP-glucose pyrophosphorylase and thereby influence the synthesis of glycogen and starch [4,11,12,14]. Therefore, research on the function of ADP-glucose pyrophosphorylase (Bna.A04g12270D) is important for improving our understanding of how plants produce and store lipids, as well as for developing new varieties of canola that have improved lipid profiles for use in food, biofuels, and other industrial applications. ACS6, Mfrn, CHUP1, LHP1, TAAT and EMB8 is not directly involved in chloroplastic fat biosynthesis, characterizing these genes is an important step in understanding the molecular mechanisms underlying fatty acid biosynthesis. Once their function is determined, these proteins can be targeted for the development of novel therapies for various diseases that are associated with dysregulated lipid metabolism.

5. Conclusions

In this study, we investigated four agronomic traits, eight fatty acid compositions, and crude fat content of the original rapeseed cultivar ‘Tamra’ and 95 rapeseed mutant lines, and found significant differences among each trait. Based on these results, upon performing a comprehensive GWAS using 73,226 filtered SNPs derived from GBS, we were able to identify 32 candidate genes significantly associated with the thousand seed weight, as well as three fatty acid compositions (C16:1, C20:1, and C22:1) and crude fat content. To strengthen the genetic basis of thousand seed weight, fatty acid composition, and crude fat content in rapeseed, future research could entail functional validation of the candidate genes that were identified in this study. Our results can offer valuable insights into the genetic basis of these traits and have the potential to info future rapeseed breeding strategies.