Genome-Wide Association Study of Fruit Traits Using 109 Germplasm Accessions of Camellia oleifera

1. Introduction

Camellia oleifera Abel, native to China, is a species within the genus Camellia of the family Theaceae. It is primarily a shrub but can also grow into a medium-sized tree. The seeds of C. oleifera possess rich contents of edible oil, also known as tea oil, making it one of four major edible oil producing crops worldwide [1,2,3]. Tea oil is rich in monounsaturated and polyunsaturated fatty acids, with significant uses in regulating blood cholesterol. Moreover, the oil contains significant contents of bioactive substances such as vitamin E and flavonoids, making it highly nutritious [4,5,6]. With a growth cycle spanning over a hundred years, C. oleifera represents a highly successful perennial cash crop.

The plant shows abundant variations in fruit shape and color among different germplasms. The most prevalent fruit colors include red, yellow, and green, while the fruit shapes can be categorized as spherical, orange-like, peach-like, navel-like, olive-shaped, jujube-shaped, and chicken-heart-shaped [7]. The key fruit traits affecting the oil production include fruit weight, fruit height, fruit diameter, fruit peel thickness, fresh seed yield rate, dry seed yield rate, and dry kernel yield rate [8,9]. These fruit traits vary among varieties and are affected by environment variations [10,11,12,13]. The diversity of fruit traits serves as an important foundation for selective breeding. In-depth analysis of the genetic basis underlying the diversity of fruit traits in C. oleifera germplasm and the identification of gene loci regulating fruit traits play crucial roles in the development of new varieties and molecular-assisted breeding [14,15,16].

However, due to the long breeding cycle, complex ploidy, and large genome size of C. oleifera, research on the genes regulating its fruit traits is relatively slow [17,18,19]. Genome-wide association study (GWAS) is an effective and convenient method for studying traits in plants and animals. By performing association analysis between genotypic and phenotypic data, significant loci or candidate genes can be identified [20]. Jia et al. [21] sequenced 916 foxtail millet varieties and identified a total of 2.58 million single nucleotide polymorphisms (SNPs). They also constructed a haplotype map of the foxtail millet genome using 0.8 million high-quality SNPs. Furthermore, they conducted a genome-wide association study on 47 agronomic traits of the 916 varieties under five different environmental conditions, revealing a total of 512 loci associated with these traits. Wang et al. [22] conducted GWAS on traits related to chlorophyll content using 529 rice germplasm resources and identified 46 significantly associated loci. They selected three of these associated loci and constructed an F₂ segregating population for validation. Wang et al. [23] also performed a GWAS using 219 soybean germplasm resources and detected a total of 12 quantitative trait loci (QTLs) related to net photosynthetic rate, photosynthetic rate, and stomatal conductance, which were distributed across 12 chromosomes. Zhang et al. [24] conducted GWAS on 11 agronomic traits of 312 sand pear germplasm resources from China, South Korea, and Japan, and identified 37 loci associated with 8 fruit quality traits and 5 loci associated with 3 fruit phenological traits. Among these, a candidate gene, PbrSTONE, was functionally validated to be involved in the regulation of stone cell formation [25]. These studies provide a framework for exploring the regulatory mechanisms of fruit traits in C. oleifera.

C. oleifera harbors rich genetic diversity in its germplasm resources [26]. This diversity provides a crucial gene pool for the genetic improvement of fruit traits in C. oleifera [27]. The development of high-throughput sequencing technology and the emergence of the C. oleifera genome offer an opportunity to develop SNP molecular markers covering the entire genome of C. oleifera [28]. However, compared with other fruit trees and crops, research on the regulatory mechanisms of fruit traits in C. oleifera remains limited. This study conducts a genome-wide association analysis on C. oleifera fruits for the first time and explores the genetic structure and diversity of C. oleifera populations. The objectives are to identify genetic loci related to the regulation of fruit traits and reveal potential molecular markers, providing a solid theoretical basis for the genetic improvement and molecular breeding of C. oleifera.

2. Materials and Methods

2.1. Plant Materials and Phenotypic Evaluation

The 109 C. oleifera germplasm accessions were planted in the experimental field of National C. oleifera Germplasm Resource Bank at Shaxian Shuinan State-owned Forestry Farm from Fujian Province (Supplementary File S1). These accessions originated from four provinces: Fujian (78 accessions), Jiangxi (17 accessions), Hunan (8 accessions), and Guangxi (6 accessions), with each accession assigned a unique identifier.

The weight of each C. oleifera fruit was measured using an electronic balance, and the length and diameter of the fruits were accurately measured using a vernier caliper. The peel thickness of each fruit was measured at the top, middle, and bottom with a vernier caliper, and the average value was calculated. The fruit shape index (FSI) was defined as the ratio of fruit length to fruit diameter. Microsoft Excel 2021 was used to calculate FSI and other statistical data.

After the C. oleifera fruits were peeled, the number of seeds in each fruit was recorded. The weight of the shelled seeds was measured, and the average value was calculated to determine the fresh seed yield rate. Subsequently, the shelled seeds were dried in an oven and then weighed again to calculate the dry seed yield rate. The dried seeds were then shelled, and the weight of the kernels was measured to calculate the dry kernel yield rate. Based on the FSI, different categories of C. oleifera fruit shapes were distinguished, identified, and recorded.

2.2. DNA Preparation and Sequencing

Genomic DNA was extracted from the leaves using a modified cetyltrimethylammonium bromide (CTAB) method [29]. After quantification by 1% agarose gel electrophoresis, the working DNA solutions were diluted to 100 ng/µL and stored at -20°C. Subsequently, we constructed sequencing libraries for the 109 DNA samples using the Hyper-seq method [30]. When the monoclonal detection results met the required standards, we performed 150 bp paired-end sequencing on each constructed sequencing library using the HiSeq 2500 platform.

2.3. SNP Calling and Annotation

The paired-end raw data from the simplified genome sequencing were filtered using the software fastp (version 0.20.0) [31], with the parameters set as: --adapter_sequence R1_adapter and --adapter_sequence_r2 R2_adapter, to obtain paired-end clean data for each sample. To facilitate rapid searching and positioning during subsequent variant detection, the index of the C. oleifera genome was constructed using three software tools: SAMtools [32], BWA [33], and Picard (http://broadinstitute.github.io/picard/). The data were aligned to the genome using the mem algorithm of BWA, and the results were then sorted using the sort function of SAMtools to obtain BAM files of various qualities. Subsequently, the alignment rate was calculated using SAMtools. PCR duplicates were removed from the BAM files using the MarkDuplicates tool of Picard, ultimately generating BAM files of various qualities containing the aligned results. Next, variant detection was performed on the BAM files of various qualities by chromosome using the HaplotypeCaller module of the Genome Analysis Toolkit (GATK) [34], generating corresponding gVCF files.

In the second step, the gVCF files of different germplasm resources were merged by chromosome using the CombineGVCFs module of GATK. Then, population variant detection was performed on the merged gVCF files using the GenotypeGVCFs module of GATK, obtaining VCF files for each chromosome. Finally, the VCF files for specific chromosomes were merged into a single genome-wide VCF file using the MergeVcfs tool of GATK.

In the third step, single nucleotide polymorphisms (SNPs) and insertions/deletions (INDELs) were extracted using the SelectVariants tool of GATK, and then filtered and flagged using the VariantFiltration tool of GATK (SNP filtering criteria: QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0; INDEL filtering criteria: QD < 2.0 || FS > 200.0 || SOR > 10.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0). Finally, the filtered SNPs and INDELs were extracted again using the SelectVariants tool of GATK, obtaining the final VCF files for SNPs and INDELs.

The identified SNPs were annotated using the ANNOVAR [35] software. First, a database was constructed using the FASTA file and GFF3 annotation file of the C. oleifera genome. Then, the obtained VCF files for SNPs were annotated.

2.4. Genetic structure Analysis

Prior to population structure analysis, the SNP data were filtered using the Plink [36] software based on criteria of missing rate (Geno 0.1), minor allele frequency (MAF 0.01), and linkage disequilibrium (LD) filtering (using the parameter --indep - pairwise 50 10 0.2). The filtered VCF files were then analyzed. The VCF files were sorted and converted into PHYLIP format using the SortGenotypeFilePlugin module of TASSEL (https://tassel.bitbucket.io/). Subsequently, a phylogenetic tree was constructed using the FastTree software [37] with the approximate maximum likelihood method and the generalized time - reversible (GTR) model. Principal component analysis (PCA) was performed on the experimental C. oleifera population materials based on the detected SNPs using the Plink software. Following this, the vectors for each principal component were calculated using R software [38], and the PCA scatter plot was plotted using R.

Additionally, the input file format was adjusted using Plink software to be compatible with the Admixture software [39]. The files were then uploaded into the Admixture software, and the number of subpopulations (k value) was set from 1 to 10. The most appropriate number of subpopulations was determined based on the cross-validation error (CV error) values generated during the computation, and a population genetic structure matrix was generated. This matrix consists of the genetic component coefficients for each individual in each subpopulation.

Finally, charts were created using the pophelper [40] and ggplot2 packages in R to visualize the results of the population genetic structure matrix. The nucleotide diversity index (π), population differentiation index (Fst), and Tajima's D values, which were calculated using the VCFtools software, were also visualized. Data were extracted from the computation result files, and the charts were plotted using the CMplot package in R.

2.5. GWAS Analysis

Prior to association analysis, SNP imputation was performed on the C. oleifera population using Beagle-5.5 [41]. The imputed data were subsequently filtered with Plink to retain markers with a minor allele frequency >0.05 and a missing rate <0.1, yielding a final set of high-quality SNPs for downstream analysis. Genome-wide association analysis for the 10 fruit traits was conducted across the 109 germplasm samples using the mixed linear model (MLM) in EMMAX, which involved kinship matrix calculation followed by association testing.

The genome-wide significance threshold was set at 3.56×10⁻⁸ (0.05/total number of SNPs), corresponding to -log₁₀(P)-value thresholds of 6.45 and 7.75. The results were visualized using the R package CMplot, generating both Q–Q and Manhattan plots in which the specified thresholds are indicated by blue and red lines, respectively.

2.6. Candidate Gene Screening

We defined the flanking region corresponding to significantly associated SNP markers as 100 kb (50 kb upstream and 50 kb downstream of the significantly associated SNP marker). Based on the annotation and functional annotation of SNP loci in C. oleifera, we considered the genes located within these loci as candidate genes. If a locus was simultaneously upstream and downstream of other genes, those upstream and downstream genes were also considered as candidate genes. If a locus was in an intergenic region, the nearest upstream and downstream genes were considered as candidate genes.

Candidate genes were annotated using GhostKOALA on the KEGG official website (https://www.kegg.jp/ghostkoala/) to obtain the relevant KO numbers. Then, a Python script was used for conversion to obtain pathway information at different levels. The Swiss-Prot database was downloaded from the NCBI official website (https://www.ncbi.nlm.nih.gov), and the protein sequences were compared to this database using BLAST (BLAST: Basic Local Alignment Search Tool) to obtain the corresponding GO numbers. Subsequently, further analysis was conducted using TBtools.

3. Results

3.1. SNP Variation Detection

After constructing libraries for each sample, approximately 126.7G of raw data were obtained through paired-end sequencing. Following filtering using the fastp software, approximately 122.7G of clean data were generated. Variants were called from the aligned sequence files of the 109 C. oleifera samples using GATK. After initial filtering, a total of 6,535,538 raw population SNPs were obtained, distributed across all 15 chromosomes. Chromosome 10 harbored the highest number of SNPs (549,986), accounting for 8.42% of the total, while chromosome 9 contained the fewest (304,057 SNPs, 4.65%). The SNP counts for the remaining chromosomes were as follows, in descending order: Chr01 (510,385; 7.81%), Chr12 (507,496; 7.77%), Chr07 (506,689; 7.75%), Chr03 (478,721; 7.32%), Chr02 (471,624; 7.22%), Chr05 (459,979; 7.04%), Chr11 (456,825; 6.99%), Chr13 (440,464; 6.74%), Chr04 (432,971; 6.62%), Chr14 (362,843; 5.55%), Chr06 (361,753; 5.54%), Chr15 (355,498; 5.44%), and Chr08 (336,247; 5.14%). Overall, the number of markers on each chromosome exceeded 300,000, indicating a relatively uniform distribution across the genome (Figure 1).

Further annotation was performed on the variant sites located within exon regions. The results showed that 152,181 SNPs were situated in exons. Among these, 90,876 SNPs led to non-synonymous mutations, 56,530 resulted in synonymous mutations, 4,562 were nonsense mutations, and 213 caused stop-loss mutations. These variants may have substantial effects on the traits of C. oleifera.

3.2. Genetic Structure

After filtering the missing rates, minor allele frequencies, and linkage disequilibrium (LD), 274,489 high-quality single nucleotide polymorphism (SNP) loci were retained. The population structures of the 109 C. oleifera germplasm accessions were analyzed using the ADMIXTURE software. The number of subpopulations (k) was set from 1 to 10, and the cross-validation (CV) error values were calculated for each k. The results showed that the CV error was the lowest when k=2, and it gradually increased as the value of k exceeded 2, indicating that dividing the germplasm resources into 2 groups at the genomic level was the most appropriate. The Q values for each germplasm accession were calculated for k values ranging from 2 to 10, and each accession was assigned to the subpopulation where its Q value was maximal (Figure 2A). The results indicated that when k=2, the 109 C. oleifera germplasm accessions exhibited distinct structural blocks, suggesting effective grouping. From statistical analysis, it was found that subpopulation Q1 (dark blue) contained 48 accessions, while subpopulation Q2 (light blue) contained 61 accessions. There was evidence of gene flow and mutual gene infiltration between the two subpopulations. Further statistical analysis showed that the germplasm accessions in Q1 largely coincided with those in subgroups G1, G2, and G3 of the phylogenetic tree. Similarly, Q2 corresponded to subgroups G4, G5, and G6. Therefore, the division of 109 germplasm resources into 2 groups yielded effective results.

Principal component analysis (PCA) was performed on the SNP data using Plink software. The first three principal components, PC1, PC2, and PC3, collectively explained 34.99% of the variance, with PC1 accounting for 12.69%, PC2 for 11.69%, and PC3 for 10.61% (Figure 2B). The results were visualized using R, employing the two subpopulations inferred by ADMIXTURE for grouping, The first three principal components, which were highly representative of the variation, were selected for plotting.

A phylogenetic tree for 109 C. oleifera germplasm resources were constructed using the approximate maximum likelihood method in FastTree software (Figure 2C). The results indicate that the 109 C. oleifera germplasm resources can be divided into six subgroups: G1 to G6. Specifically, subgroup G1 contains 12 germplasm accessions, G2 has 27, G3 has 17, G4 has 18, G5 has 11, and G6 has 24. In subgroup G1, 2 accessions are from Fujian, 4 are from Hunan, and 6 are from Jiangxi. In subgroup G2, 14 accessions are from Fujian, 5 belong to Guangxi, 3 are from Hunan, and 5 are from Jiangxi. Subgroup G3 includes 11 accessions from Fujian and 6 from Jiangxi. Subgroups G4 and G5 contain 18 and 11 accessions, respectively, all from Fujian. In subgroup G6, apart from 1 accession from Hunan and 1 from Jiangxi, the remaining 22 are from Fujian. Further analysis reveals that the germplasm from Fujian in subgroups G4, G5, and G6 includes 'Min 43', 'Min 46', 'Min 49', and 'Min 50', as well as their hybrid offspring 'Min Zayou'. These results suggest that the distinct geographical distribution characteristics of subgroups G4, G5, and G6 may be influenced by hybridization, leading to complex genetic backgrounds within these subgroups and subsequent population differentiation. Overall, there is no obvious geographical distribution pattern among the different subgroups of the 109 C. oleifera germplasm resources.

Further calculations were performed to determine the distribution of Fst for each SNP marker across the whole genome (Figure 3). The Fst values for most SNP markers were below 0.2, and the average Fst value for all SNP markers across the genome was 0.0153. This suggests that the degree of differentiation between the Q1 and Q2 subpopulations is relatively weak, indicating weak genetic differentiation.

To assess the genetic diversity of the two subpopulations within the C. oleifera germplasms, the nucleotide diversity index (π) was calculated for the Q1 and Q2 subpopulations using VCFtools software. The distribution of π values across the 15 chromosomes for both subpopulations is shown in Figure 4A,B. There was no significant difference in genetic diversity between the two subpopulations, and both exhibited higher nucleotide diversity at certain SNP loci on chromosomes 6 and 8, indicating that these genetic variations have been well-preserved.

In addition, the average nucleotide diversity index for the Q1 subpopulation was 4.45×10-4. The nucleotide diversity indices for both subpopulations were mostly below 0.01, depicting low genetic diversity. This may be due to long-term domestication of C. oleifera.

Finally, the Tajima's D values were calculated for the two subpopulations (Figure 4C,D). Tajima's D is an indicator of selection within populations and is used to identify sequences that deviate from the neutral theory model by assessing the balance between mutation and genetic drift. The mean Tajima's D value for the Q1 subpopulation across all SNP loci was 0.292, while it is 0.464 for the Q2 population. These positive Tajima's D values suggested that the presence of rare alleles at low frequencies in both subpopulations, and there was a relatively high frequency of intermediate alleles. This may be due to the influence of domestication on the C. oleifera germplasm, leading to similar environments and uniform selection pressures within the populations.

3.3. GWAS of Fruit Traits

A total of 2,812,326 SNP markers spanning the C. oleifera genome and 10 fruit traits (Supplementary File S2) from 109 germplasms were utilized for GWAS analysis. The analysis was performed using the EMMAX mixed linear model, with significance thresholds set at 6.45 and 7.75 for screening significant loci. All traits except for fruit height were found to be associated with significant loci. In total, 157 significant loci were identified (Figure 5).

Specifically, nine SNP loci were significantly associated with single fruit weight, among which the locus at Chr09:7869152 showed a highly significant association (-log₁₀(P) = 8.05). Two loci associated with fruit diameter were located at positions 43866733 and 43866738 on chromosome 5. For the fruit shape index, six significant loci were detected, with the most significant at Chr02:20212286 (-log₁₀(P) = 7.42). Forty loci were associated with peel thickness; 37 showed significant and three showed highly significant associations. The SNP at Chr13:111073692 was the most significantly associated (-log₁₀(P) = 10.01), followed by Chr10:33172723 (8.17) and Chr11:149010807 (7.76). Among the ten traits, the number of seeds per fruit had the highest number of associated loci (69 in total: 59 significant, 10 highly significant). The top three loci were Chr11:153463775 (-log₁₀(P) = 10.09), Chr15:29597100 (9.42), and Chr15:29597154 (9.01). Five loci were associated with fresh seed weight, with the most significant at position 169280653 on chromosome 10 (-log₁₀(P) = 6.91). Twenty-three loci were associated with the fresh seed yield rate, with only one (Chr11:87995439) being highly significant (-log₁₀(P) = 8.56). Only one locus (Chr07:160272823) was associated with the dry seed yield rate (-log₁₀(P) = 6.68). Finally, two loci (Chr02:124313904 and Chr02:124313920) were associated with the dry kernel yield rate, both with a -log₁₀(P) value of 6.72 (Supplementary File S3).

3.4. Candidate Genes

Annotations were performed on the 157 significant SNP loci to identify for candidate genes. After removing duplicate genes, 110 candidate genes (Supplementary File S3) were annotated using the KEGG database (Table 1). The KEGG annotations primarily included genetic information processing protein families, metabolic protein families, terpenoid and polyketide metabolism, among others.

The GO annotations of 46 genes (Table 2) suggested that 6 candidate genes belonged to the molecular function category, including other molecular functions such as disease resistance. Nine candidate genes were associated with the cellular component category, mainly involving functions such as organelle membranes, protein complexes, and the endoplasmic reticulum. Thirty-one candidate genes showed biological processes, primarily involving functions such as reproduction, development, and RNA biosynthesis.

4. Discussion

In recent years, SNP molecular markers have assumed a critical role in elucidating biological genetic variation, interspecies relationships, and genetic breeding, owing to their numerous advantages in comparison to traditional breeding methods. These advantages include their abundance, wide distribution, ease of automated detection, low mutation rate, and rapid screening [42,43]. Currently, Simple Sequence Repeat (SSR) [44], Inter Simple Sequence Repeat (ISSR) [45], and Sequence Related Amplified Polymorphism (SRAP) [46] molecular markers are commonly used to study the genetic diversity of C. oleifera. However, these markers have limitations, such as high time costs, low throughput, and relatively low accuracy, making it difficult to screen for highly polymorphic markers. Liu et al. [47] utilized SLAF-seq technology to sequence nine C. oleifera samples and successfully developed SNP markers. However, we did not use the C. oleifera genome as a reference for variant detection; consequently, marker development requires improvement and is insufficient for gene mapping of C. oleifera traits. Furthermore, the population size used in the study was relatively small, rendering the quality of the developed markers inadequate for comprehensive investigation of genetic diversity within C. oleifera populations.

With the rapid development of bioinformatics and continuous advancements in sequencing technology, the genome of C. oleifera has been successfully sequenced. Lin et al. [48] conducted transcriptome sequencing using the C. oleifera genome as a reference and detected 1,849,953 SNP molecular marker loci. Based on these SNP loci, they constructed high-quality molecular identities for 221 common C. oleifera resources. In this study, we performed ddRAD reduced-representation genome sequencing on 109 C. oleifera germplasms, using the C. oleifera genome as a reference and GATK4 for variant detection. A total of 6,252,197 population SNP loci and 722,237 population INDEL loci were identified. Compared to the study by Lin et al. [48], our study had a higher density of SNP molecular markers. Using ANNOVAR software for annotation analysis of the variant loci, we found that 145,679 SNPs and 12,756 INDELs were in exon regions. Variants in exon regions may lead to changes in the amino acids encoded by codons [49], and even mutations in start and stop codons, thereby altering the function of polypeptide chains and proteins and further affecting the traits of biological individuals. Therefore, variant loci located in exon regions have a significant impact on biological traits, and these variants may significantly influence the traits of C. oleifera.

Based on the combined results of Structure analysis, phylogenetic analysis, and principal component analysis, it is reasonable to divide the 109 C. oleifera germplasms into two groups: Q1 and Q2. The average Fst value across all SNP markers in the entire genome between these two groups is 0.0153, indicating a very weak level of differentiation between the two populations. Lin et al. [48] classified 221 C. oleifera materials into eight subgroups based on genetic structure, with Fst values ranging from 0 to 0.0085 among the subgroups. These findings further indicate a weak level of differentiation among the C. oleifera subgroups. Our results are consistent with their findings. Furthermore, Huang [50] used SSR molecular markers to analyze C. oleifera populations and found that the genetic differentiation coefficient of populations is lower than that of woody angiosperms. Therefore, the weak differentiation of C. oleifera populations may be attributed to the fact that it is a cross-pollinated plant that has undergone long-term domestication and extensive hybridization, leading to high gene flow among subgroups.

The average nucleotide diversity index for the Q1 subgroup is 4.45×10-4. These results are generally consistent with those of Lin et al. [48]. However, compared to Brazilian cassava plants (π=0.274) [51], the nucleotide diversity index of C. oleifera populations is very low, indicating a relatively low level of genetic diversity in germplasms. This may be due to long-term domestication of C. oleifera [52]. To obtain high-yielding and high-quality C. oleifera varieties, plants have been subjected to long-term domestication, targeted breeding, and continuous selection, resulting in improved agronomic traits. However, long-term domestication has favored alleles that establish population advantages, while relatively disadvantageous alleles have been eliminated. Over time, this has led to a decrease in genetic diversity within C. oleifera populations [48]. The Tajima's D values for both subgroups are greater than zero, indicating the presence of rare alleles at low frequencies and a relatively high number of intermediate-frequency alleles in both subgroups. This may be due to the similar growth environments of C. oleifera, leading to uniform selection pressure across the population.

The whole genome association analysis (GWAS) is based on linkage disequilibrium and involves correlation analysis between target traits and variation sites to identify significant variation sites, thereby uncovering candidate genes for the target traits. GWAS has made significant achievements in the mapping of complex traits and the improvement of crop varieties in both plants and animals, such as rice [53], apple [54], and walnut [55]. However, due to the gaps in the sequencing and annotation of the C. oleifera genome, there have been few studies using GWAS to analyze C. oleifera traits. We identified 157 genetic loci significantly associated with fruit traits in C. oleifera. Specifically, there were 9 significant loci for single fruit weight, 2 for fruit diameter, 6 for fruit shape index, 40 for peel thickness, 69 for seed number per fruit, 5 for fresh seed weight, 23 for fresh seed yield rate, 1 for dry seed yield rate, and 2 for dry kernel yield rate. Notably, no significant loci were found for fruit height. We further screened 110 candidate genes based on the significant SNP loci. Among these, 46 candidate genes showed GO annotations. Among them, 6 candidate genes belonged to molecular function, 9 to cellular component, and 31 to biological process, involving various gene functions. Remarkably, 16 out of the 31 candidate genes related to biological processes were involved in reproduction, indicating their important roles in the development of C. oleifera fruits. In the Kyoto Encyclopedia of Genes and Genomes (KEGG) annotation, 23 genes were successfully annotated, primarily associated with genetic information processing protein families, metabolic protein families, and terpenoid and polyketide metabolism. These significant SNP loci and candidate genes may be involved in multiple traits such as yield, fruit size, and peel thickness, providing important genetic information and molecular markers for the improvement of quality and yield of C. oleifera.

In our study, the annotation results of significant SNP sites indicate that most SNP markers are in intergenic regions. This may be attributed to the widespread distribution of intergenic regions in the C. oleifera genome [56]. Intergenic regions were once considered non-functional; however, the completion of sequencing for an increasing number of plant genomes, it is revealed that intergenic regions contain open chromatin. Moreover, it has been observed that larger plant genomes exhibit a higher proportion of open chromatin within their intergenic regions. Chromatin accessibility is involved in plant environmental responses, growth, and development [57,58]. The functions of intergenic regions in the C. oleifera genome require further research and exploration. Therefore, candidate genes located both within genes and in intergenic regions that harbor significant SNP sites have the potential to influence corresponding traits.

This study offers novel insights and empirical data to support research on the genetic evolution and genomic structure of C. oleifera. Moving forward, the integration of significant SNP loci and candidate genes will facilitate more comprehensive functional studies and gene editing experiments, thereby elucidating the molecular mechanisms underlying key traits of C. oleifera. This will provide more reliable and effective technical support for the genetic breeding and industrial development of C. oleifera.

5. Conclusions

This study performs a genetic diversity analysis of 109 C. oleifera germplasm accessions using fruit phenotypic traits and SNP molecular markers. GWAS identified significant SNP loci and candidate genes influencing fruit phenotypic traits. By integrating structure analysis, phylogenetic analysis, and PCA, we classified the 109 C. oleifera germplasm accessions into two distinct populations. The genetic differentiation between the two subgroups was weak, with nucleotide diversity indices (π) of 4.45×10-4 and 5.35×10-4, respectively, indicating low genetic diversity in C. oleifera germplasms. This may be attributed to long-term domestication of the species. The Tajima's D values for both subpopulations were greater than zero, suggesting the presence of numerous intermediate-frequency alleles. This may be attributed to the similar environments in which the C. oleifera germplasm grows, leading to uniform selective pressures within the population. GWAS was employed to pinpoint loci significantly associated with fruit traits in the germplasm resources. The proteins encoded by the identified candidate genes were primarily involved in various biological processes, including biosynthesis, regulation of stem system morphogenesis, reproduction, plant organ development, transport, and metabolic pathways. These proteins include plant hormone glycosyltransferases, peptidases and inhibitors, ubiquitin protein ligases, cytochrome P450s, and protein kinases. Notably, among the candidate genes related to 31 biological processes, 16 were specifically involved in reproduction. This highlights the crucial role these candidate genes play in the development of C. oleifera fruits. Our findings provide novel insights and data support for the study of genetic evolution and genomic structure in C. oleifera. Furthermore, they hold significant importance for the targeted breeding of elite varieties, paving the way for future improvements in C. oleifera cultivation.