Association Mapping of Hundred-Grain Weight in Paeonia ostii Using SSR Markers of Transcription Factors from the Comparative Transcriptome

1 Introduction

Tree peony (Paeonia sect. Moutan DC), a perennial deciduous shrub native to China, is highly valued for its significant horticultural and medicinal properties. It is widely acknowledged as both a prestigious ornamental plant and an indispensable component of traditional Chinese medicine. In 2011, China’s Ministry of Health recognized tree peony seed oil as a new food resource. Subsequently, in 2014, the China Food and Drug Administration further incorporated it into the list of approved raw cosmetic materials, thereby officially granting tree peony seed oil enter onto the market for edible oil and cosmetics [1].

P. ostii is a critical germplasm resource for oil peony cultivation. The seed oil of P. ostii contains more than 90% unsaturated fatty acids, with approximately 42.78% α-linolenic acid (ALA), approximately 26.50% linoleic acid (LA), and approximately 22.97% oleic acid (OA) [2]. Extensive research has demonstrated that ALAs offer substantial health benefits, including antiglycation effects [3], the promotion of brain development [4], the protection of cardiovascular health [5], neuroprotective functions [6], and a reduction ininflammatory responses [7].

However, since ALA is an essential fatty acid that mammals cannot synthesize de novo, humans must rely on dietary intake to meet their physiological requirements [8]. Thus, P. ostii seed oil is positioned as a nutraceutical-grade oil, combining health-promoting attributes with promising economic viability. Moreover, as one of the peony species that exhibit exceptional adaptability, P. ostii has demonstrated remarkable tolerance and resilience to diverse environmental conditions. This characteristic renders it highly suitable for cultivation in nonarable land areas, such as mountainous regions, hilly terrains, and forest margins, thereby effectively conserving farmland resources [1]. In conclusion, as an important woody oil crop, P. ostii not only addresses China’s national strategic needs for producing healthy edible oil, promoting economic development, and ensuring ecological protection but also has remarkable application value and extensive development potential. This topic merits further in-depth research and broad promotion throughout China.

To date, research on the seeds of P. ostii has focused predominantly on specific aspects, such as the characterization of fatty acid composition in seed oil [2], the evaluation of nutritional components [9], and the elucidation of fatty acid metabolic pathways [10,11]. Prior studies have indicated that the oil content of P. ostii seeds is approximately 30%, which is considerably lower than that of conventional oilseed crops such as rapeseed (35–50%), sunflower (40–50%), and peanut (46–57%) [12]. Additionally, it is lower than that of some emerging woody oil plants, including Xanthoceras sorbifolium (24–55%) [13], Camellia oleifera (41–53%) [14], and Juglans regia (54–74%) [15].

The relatively low oil content in peony seeds directly leads to a reduced yield of tree peony seed oil, which subsequently increases its market price and reduces consumer demand. Consequently, enhancing the oil content and yield of P. ostii seeds is crucial for advancing their application in edible oil production. Nevertheless, substantial variations exist in the oil content and yield of P. ostii seeds, both across diverse geographical regions and among individuals within the same region [16,17]. The reproductive strategy of tree peony relies primarily on outcrossing, with its cultivated diversity largely driven by hybridization [18,19].

The genetic improvement of tree peony germplasm traditionally depends on controlled crosses and phenotypic selection. The development of a new variety with stable traits typically requires at least 10 years. Accelerating genetic improvement and seed output through molecular-assisted selection is now imperative. To expedite the development of high-yield varieties, it is essential to conduct an in-depth investigation into the molecular mechanisms underlying tree peony seed yield formation. The initial phase of this research could center on analyzing the critical factors influencing peony seed yield. Crop yield, a multifaceted quantitative trait governed by environmental and genetic interactions, depends primarily on three key determinants: grain weight, panicle count, and grains per panicle [20,21,22,23]. Among the three elements that constitute crop yield, compared with the panicle count and number of grains per panicle, the grain weight has is more heritable and has more stable genetic characteristics. Therefore, in research on improving crop yield, grain weight is regarded as a more crucial and influential trait [24,25].

Grain weight is primarily determined by grain size, which encompasses grain length, width, thickness, and length-to-width ratio (LWR). Key factors influencing grain size include the ubiquitin-proteasome pathway, G-protein signaling, the mitogen-activated protein kinase (MAPK) signaling cascade, plant hormone regulatory mechanisms, and transcriptional regulator activities. Moreover, the efficiency of nutrient utilization and the regulation of storage product accumulation play crucial roles in determining grain weight. These mechanisms collectively influence the final grain weight [26]. However, currently, research reports concerning the grain weight of tree peony remain relatively limited, presenting significant challenges for the direct selection of peony seed weight through molecular approaches. Consequently, analyzing the genetic basis of yield traits on the basis of a molecular marker linkage map, finely mapping the quantitative trait loci (QTLs) that control grain weight, and thoroughly investigating the molecular genetic mechanisms underlying tree peony seed weight formation are crucial steps toward accelerating the high-yield breeding process of P. ostii.

The inheritance of quantitative traits typically involves multiple genes, each exerting a relatively small effect and being susceptible to environmental influences. These traits are generally characterized by low heritability and have historically been considered challenging to study. Nevertheless, the rapid advancement of molecular marker technologies, genomic mapping approaches, and QTL analysis has facilitated significant progress in elucidating the genetic architecture underlying quantitative traits. Owing to the vastness and complexity of the tree peony genome [27], current population genetics analyses of tree peony still rely mainly on microsatellite markers (SSRs). In plant breeding, microsatellite markers serve as crucial tools for marker-assisted selection (MAS) and continue to play an indispensable role. Their advantages, including codominant inheritance patterns, high transferability across species, strong associations with functional genes, and relatively low development costs, have made SSR markers highly attractive in the field [28]. In particular, coding-region SSR loci, which function as expressed transcription factors, demonstrate significant evolutionary conservation across species and directly modulate gene transcription and translation processes [29]. The substantial accumulation of transcriptome data in tree peony provides a solid foundation for the development of SSR markers [30,31,32,33,34,35,36,37,38].

To date, more than 500 SSR markers have been developed for tree peony [39,40,41,42,43,44,45], and these markers have been widely used in research fields such as genetic diversity assessment, population structure analysis, cultivar identification, and genetic map construction. However, the coverage of existing molecular markers, especially EST-SSR markers, remains relatively limited. This situation has, to a certain extent, constrained broader applications of functional gene markers across diversity assessment, genetic clustering characterization, phylogenetic population segmentation, linkage disequilibrium research, and association analysis.

In conclusion, research on functional genes associated with yield traits in P. ostii remains relatively limited, particularly regarding the genetic architecture of grain weight—a key yield component. To date, no association mapping has been specifically conducted for hundred-grain weight in this emerging woody oil crop. In this study, we leveraged this resource by integrating long-read and short-read sequencing to perform the first mapping association for hundred-grain weight in P. ostii. First, on the basis of EST-SSR markers developed through comparative transcriptomics, we systematically evaluated the genetic diversity and population structure. Subsequently, by performing linkage disequilibrium analysis and single-marker association mapping, we further screened the effects of key alleles on the natural variation in grain weight. Our findings not only reveal novel markers and candidate genes for marker-assisted selection but also provide a foundation for understanding the molecular mechanisms of seed weight regulation in woody perennials.

2 Methods

2.1 Plant Materials

In accordance with the principle of maximizing the coverage of existing phenotypic variation, a total of 123 individuals were sampled across most of the distribution area in China. These samples originated from Henan Province in central China, a primary cultivation region for P. ostii, and were cultivated in the same nursery field at Henan Agricultural University in Zhengzhou, China (34°51′ N, 113°35′ E), for more than seven years under randomized conditions. All the evaluated plants exhibited normal flowering and fruiting patterns and presented stable and consistent traits. Overall, the plants showed more consistent year-to-year performance in terms of both flowering and fruiting.

2.2 Phenotypic Data

In September 2024, at the full maturity stage (characterized by pericarp dehiscence), the seeds were harvested and air-dried uniformly in a controlled environment (25°C, 30% relative humidity) for 72 h to a constant moisture content. The hundred-grain weight was measured using a precision electronic balance (METTLER TOLEDO MS104TS, readability: 0.001 g) in a temperature-stable laboratory. Three biological replicates were measured per accession. The mean value derived from three replicates per individual was utilized for statistical and association analyses. Descriptive statistics for hundred-grain weight, including the minimum, maximum, mean, and standard deviation (SD), were calculated using SPSS 18.0 (IBM Inc., Chicago, IL, USA). The coefficient of variation (CV) was subsequently derived as the (SD/mean) × 100%. Individual variability was assessed through one-way ANOVA. Prior to analysis, all SM datasets were subjected to normality verification via the Kolmogorov-Smirnov (K-S) test and variance homogeneity evaluation using Levene’s test. Statistical significance thresholds were set at p < 0.05 (significant) and p < 0.01 (highly significant) for differential interpretation.

2.3 EST-SSR Selection

We performed comparative transcriptomic profiling of carpel quantitative traits in two P. rockii cultivars (‘Jing shun fen’ vs. ‘Fen mian tao sai’) through combined long-read sequencing (PacBio) and short-read (Illumina) RNA-seq data [46]. Sequencing analysis revealed 959 candidate transcription factors across 21 yield-associated gene families, alongside 166 EST-SSRs encompassing six nucleotide repeat motifs. Among these, 57 polymorphic EST-SSR markers (Table S1) were validated as effective in Paeonia [47]. In addition to the 24 previously reported polymorphic EST-SSR markers (Table S2) related to tree peony seed yield [48], 81 primer pairs were employed for genotype screening across 123 accessions.

2.4 DNA Isolation with EST-SSR Genotyping

Genomic DNA was isolated from 30 mg silica-desiccated foliage using a DNAsecure Plant Kit (Tiangen Biotech, Beijing, China). Purity and quantification were determined through 260 nm absorbance measurements with an Agilent UV–Vis spectrophotometer. Stock solutions were standardized to 30–50 ng/μL and cryopreserved at −20°C for downstream applications.

SSR–PCR amplification was performed in 10 μL reactions containing 5 μL of 2× Power Taq PCR Master Mix (Aidlab Biotechnologies, Beijing, China), 3 μL of distilled water, 1 μL of genomic DNA (30–50 ng/μL), and 0.5 μL of each forward and reverse primer (10 pmol), which were labeled with fluorescent dyes such as 6-FAM, HEX, TAMRA, or ROX (Ruiboingke, Beijing, China).

PCR amplification proceeded as follows: 5 min initial denaturation (95°C); 35 cycles of 95°C/30 s denaturation, Ta-specific annealing (48.7–61°C; Supplementary Tables S1 and S2)/30 s, and 72°C/1 min extension; concluding with a 10 min final extension at 72°C.

Amplified fragments were subjected to capillary electrophoresis on an ABI3730XL system (Applied Biosystems, Carlsbad, CA, USA) with LIZ 600 standards. Polymorphic EST-SSR alleles were resolved through GeneMarker v2.2 (SoftGenetics, State College, PA, USA) analysis.

2.5 Population Genetic Diversity and Structure

Population genetics parameters (NA, NE, I, HO, HE, FIS, and GD) were calculated using GenAIEx version 6.5, while the Microsatellite Toolkit was used to assess marker-specific PIC values. Three methods were used to analyze the population structure on the basis of 81 EST-SSR markers. First, genetic clustering analysis implemented in STRUCTURE v2.3.4 inferred population subdivisions (K = 1–10) through 10 replicate runs per K-value, featuring 100,000 burn-in/200,000 sampling iterations. The program implements Bayesian model clustering through Markov chain Monte Carlo (MCMC) integration. The results of these runs are subsequently uploaded to StructureSelector [49] to determine the final optimal K and to visualize the optimal number of clusters. The optimal K was determined using both LnP(D) and ΔK [50]. The matrix corresponding to this optimal K can be utilized to adjust for false positives in single-marker analyses.

The second approach utilized neighbor-joining (NJ) clustering of Nei’s unbiased genetic distances computed through MEGA-X to construct phylogenetic trees [51]. For the third analytical framework, principal coordinate analysis (PCoA) was performed using GenAIEx v6.5, implementing the same Nei’s distance metric across all the germplasm accessions.

2.6 LD and Phenotype—Genotype Association Analysis

The level of LD between different loci was quantified using the squared allele frequency correlation coefficient (r2). An r2 threshold of 0.1 was set to identify significant linkage disequilibrium (r2 > 0.1). TASSEL v2.0.1 was used to calculate the r2 between SSR marker pairs through 105 iterative analyses. Pairs of SSR markers were deemed to be in significant linkage disequilibrium when p was less than 0.001.

Single-marker association analysis for hundred-grain weight was implemented with EST-SSR markers through the MLM algorithm in TASSEL v2.0.1. The MLM incorporates both population structure and the kinship matrix to minimize false-positive associations. In this context, Q represents the matrix derived from the optimal number of subgroups identified through STRUCTURE analysis; K denotes the kinship matrix reflecting relationships among individuals, which is calculated using SPAGeDi v1.2. Subsequently, association mapping was conducted by combining the phenotypic database, genotypic database, Q matrix and K matrix data. Finally, the false discovery rate (FDR) for the p-value (p < 0.01) of the association mapping results was evaluated using the R package QVALUE. Markers with an FDR-adjusted Q-value < 0.05 were considered significantly associated with the hundred-grain weight.

Significant loci underwent gene action characterization through dominance-to-additivity (d/a) coefficient ratios. Classification thresholds were established as follows: |d/a| ≤ 0.50 (additive), 0.50 < |d/a| < 1.25 (partial/complete dominance), and |d/a| > 1.25 (overdominance). The algorithms and formulas used for these calculations were detailed in previous studies [52].

Full-length sequences of candidate transcription factor genes were obtained from the P. ostii transcriptome assembled from combined PacBio long-read and Illumina short-read sequencing data. Complete open reading frame (ORF) prediction and translation were performed on candidate markers using Open Reading Frame Finder (ORF Finder), and sequence identities were confirmed by BLASTx against the NCBI nonredundant protein database. The resulting proteins were screened against the Arabidopsis Information Resource (TAIR) for homology analysis, followed by phylogenetic reconstruction and sequence alignment visualization through DNAMAN.

3 Results

3.1 Genetic Diversity

One-way ANOVA of the hundred-grain weight revealed highly significant differences among individuals (p < 0.01) (Table S3), with weights ranging from 20.9 to 77.07 g and a coefficient of variation of 20.33%. Statistical validation demonstrated homogeneity of variance (F-test, p = 0.67) and normally distributed data (Kolmogorov-Smirnov, p = 0.20). Furthermore, both the distribution histogram (Fig. S1) and the probability-probability (P-P) plot (Fig. S2) for the hundred-grain weight further illustrate that the data were normally distributed, fulfilling the prerequisites for further analysis.

Genetic diversity analysis of 123 germplasm accessions using 81 polymorphic EST-SSR markers revealed 467 alleles (2–14 per locus; mean 6.0). Key parameters included effective alleles (1.008–7.952; avg 2.585), Shannon’s index (0.026–2.260; avg 1.018), heterozygosity (HO: 0.008–1.000; avg 0.501; HE: 0.008–0.874; avg 0.522), and FIS coefficients (−0.953–0.821; avg 0.040) for 46 loci, with HO > HE and PIC values (0.008–0.862; avg 0.469) indicating species-level diversity (Table 1).

3.2 Population Structure

Population structure, which refers to the presence of subgroups within a population, is a critical factor influencing LD. The existence of these subgroups can intensify LD, leading to an increased likelihood of false positive results. Consequently, it is essential to analyze and adjust for population structure prior to conducting association studies.

STRUCTURE analysis revealed that LnP(D) tended to gradually rise as K increased, but no distinct peak was observed (Fig. 1a). We subsequently statistically analyzed the range of ΔK for K = 1 to 10. On the basis of the principle that the optimal number of subgroups corresponds to the K with the maximum ΔK, we determined that the optimal number of subgroups K for P. ostii in this study was 2 (Fig. 1b). On this basis, we conducted a detailed analysis of the Q matrix and the gene frequency profiles of individual samples for K values ranging from 2 to 8 (Fig. 1c). The results demonstrated that the Bayesian method identified extensive numbers of heterozygous lineages within the population. When K ranged from 2 to 7, the population samples were consistently divided into two primary clusters. Specifically, starting from K = 2, the population samples exhibited two relatively stable clusters, suggesting that the optimal subgroup structure of the population was two.

Figure 1: Estimation of genetic structure of 123 accessions for P. ostii population using 81 SSRs based on the STRUCTURE. (a) Log probability data [LnP(D)] for each K value (10 replicates). (b) ΔK estimates of the posterior probability distribution of the data for a given K. (c) Estimated population structure and clustering of the 123 P. ostii individuals with K = 2 to 8. Individuals are shown by thin vertical lines, which are divided into two major well-separated genetic clusters standing for the estimated membership probabilities.

On the basis of the allele frequencies, genetic distances among individuals within the population were calculated. Using the NJ method, a phylogenetic tree for the associated population was constructed. These individuals initially clustered into two major clades (Fig. 2), which closely mirrored the classification results obtained from the STRUCTURE analysis. As a complementary approach to STRUCTURE-based clustering, principal coordinate analysis with SSR-derived eigenvectors revealed two genetically distinct subgroups within the cultivated P. ostii population (Fig. S3). This stratification pattern was consistent across the three independent analytical frameworks.

Figure 2: The phylogenetic tree of association population based on 81 SSRs.

3.3 LD Level

The level of LD was evaluated for 81 SSRs across 123 individuals (Fig. 3). The squared correlation coefficient of allele frequencies (r2) was used to measure the extent of LD between markers. The analysis revealed that 81 EST-SSR markers produced a total of 3240 pairwise comparisons.

Figure 3: Pairwise LD (r2) between SSRs. X and Y axis represent the 81 SSRs. The different colors correspond to the thresholds of r2 and p. r2 < 0.1 represent linkage equilibrium, r2 > 0.1 represent linkage disequilibrium and p represents the significant difference.

The r2 values ranged from 0.00006719 to 0.62929244. Notably, only 0.31% (r2 ≥ 0.1) of these comparisons showed significant LD, demonstrating predominant linkage equilibrium across loci (r2 < 0.1; p < 0.001). For instance, within the MYB gene family, three markers—PS10, PS12, and PS17—were in linkage equilibrium. However, several markers displayed significant LD (r2 > 0.5; p < 0.001), such as PS50, within the gene annotated as PsTCP11.

3.4 Single-Marker Associations of Hundred-Grain Weight in P. ostii

The MLM revealed 81 marker-trait associations for hundred-grain weight, including 21 significant hits (25.93%) at p < 0.01 (Table 2). The single marker explained 9.15% to 35.94% (with an average of 18.61%) of the phenotypic variation. FDR correction based on multiple testing was subsequently applied to 21 significantly associated combinations, and it was found that all pairs remained significantly associated (Q < 0.05). Gene effect quantification revealed additive (6, 28.57%; |d/a| < 0.5), dominant (2, 9.52%; 0.5 < |d/a| < 1.25), and over-dominant (13, 61.90%; |d/a| > 1.25) modes among the 21 trait-marker associations.

3.5 TFs Associated with Hundred-Grain Weight in P. ostii

Association analysis revealed that the TFs significantly associated with hundred-grain weight belonged to 11 gene families (AP2, TCP, YAB, WRKY, WD, BHLH, MADS-box, etc.) (Table 2). Furthermore, among the significantly associated genes, two genes were identified from the AP2, TCP, and YAB gene families. Additionally, the genes with an explanation rate exceeding 20% encompassed PS147 (WD), PS95–96 (YAB), PS163 (TFIID), PS27 (AP2), P280 (AP2), and PS98 (B3). Through sequence alignment via the National Center for Biotechnology Information (NCBI) and the Arabidopsis Information Resource (TAIR), the gene sequences of P280, PS27, PS95–96, and PS98 were found to be associated with growth and development.

P280 was predicted to encode WRINKLED1 (WRI1), a transcription factor characterized by a typical APETALA2 (AP2) domain, which plays a crucial role in the regulation of storage compound biosynthesis. After sequence alignment, it was found that P280 corresponded to PoWRI1 in P. ostii. PoWRI1, which was localized in the nucleus and associated with the oil accumulation process, had been cloned and identified. Overexpression of PoWRI1 in Arabidopsis resulted in larger seeds in transgenic plants, a significant increase in oil content, and an increase in unsaturated fatty acid levels [53]. Furthermore, in Arabidopsis, AtWRI1 was found to exhibit 53% similarity to PoWRI1 (Fig. S4), and its primary function was also widely recognized as regulating oil biosynthesis [54]. PS27 was predicted to be AtERF7, encoding an ERF subfamily B-1 protein (ethylene response factor) within the AP2/ERF transcription factor family (Fig. S5). In Arabidopsis, WRKY22, GATA8, and ERF7 collectively play critical roles in nitrogen metabolism, which in turn regulates plant growth and development [55].

The explanatory rate of PS95–96 was 33.04%, which was attributed to YAB1 (Fig. S6). YAB1, a YAB family transcription factor, phylogenetically regulates plant morphogenesis through conserved developmental pathways. More specifically, it is indispensable for meristem development and inflorescence formation [56].

The predicted result of PS98 was AtNGAL3, which encodes a plant-specific B3 DNA-binding domain transcription factor. PS98 was annotated as AtNGAL3, encoding a plant-specific transcription factor containing a B3 DNA-binding domain (Fig. S7). AtSOD7/AtNGAL2 encodes a transcription repressor belonging to the B3 family, specifically NGAL2 (NGATHA LIKE2). The overexpression of SOD7/NGAL2 leads to reduced seed and organ size, whereas the simultaneous knockout of SOD7/NGAL2 and its closest homolog, DPA4/NGAL3, results in significant increases in seed and organ size. These findings suggest that SOD7/NGAL2 and DPA4/NGAL3 redundantly regulate seed and organ size [57].

4 Discussion

4.1 Population Genetic Variance

An essential prerequisite for conducting association analysis is to precisely elucidate the genetic diversity within the associated population and the polymorphism level of the markers used [58]. Here, we evaluated the hundred-grain weight of an association population consisting of 123 individuals of P. ostii. The results revealed that the average weight of in the 123 P. ostii individuals was 45.39 g, which was close to the 48.10 g reported by Han et al. [59] and the 42.69 g/41.78 g determined by Wang et al. [60,61], and greater than the 37.14 g/37.21 g/38.36 g of P. rockii [62]. The coefficient of variation (CV) quantifies normalized dispersion through standard deviation-to-mean standardization, serving as an indicator to reflect the relative extent of data dispersion while effectively eliminating the influence of measurement scale and dimensionality. In general, a larger CV indicates a higher degree of trait value dispersion. The overall CV for each phenotypic trait is calculated on the basis of the average CV across groups, ensuring comparability between groups and avoiding interference from factors such as plant age. In this study, the hundred-grain weight of P. ostii ranged from 20.90 g to 77.07 g, with a CV of 20.33%. According to Tan et al. [63], the hundred-grain weight of 239 P. ostii individuals varies between 7.00 g and 70.00 g, with a CV of 20.30%. On the one hand, the extremely close CVs in both studies suggest that the degree of seed weight variation in P. ostii is highly consistent across different investigations. On the other hand, the significant differences among individuals in P. ostii can be attributed to its wide geographical distribution, insect-mediated pollination mechanism, outcrossing reproduction mode and self-incompatibility breeding system. These phylogenetic signatures collectively enable this species to maintain relatively high genomic variation, which is in accordance with that of other allogamous species [64].

Analysis of phenotypic variation enables a preliminary understanding of the characteristics of associated populations. However, as phenotypic variation is highly susceptible to environmental factors and the phenotypic differences among certain individuals are minimal, relying solely on phenotypic traits makes it challenging to achieve effective differentiation. Therefore, assessing genetic variability using molecular markers is particularly important. In this research, 123 P. ostii individuals were analyzed using 81 EST-SSRs, revealing 467 alleles with an NA of 5.765. This value was significantly greater than those reported in previous studies 1.900 (NA) based on 22 SSRs for 15 wild P. ostii individuals [65] and 2.259 (NA) based on 29 SSRs for 901 cultivated P. ostii individuals [66]. Compared with those in previous studies examining other species of Paeonia sect. Moutan, the mean allele count in this study was greater than that observed in P. delavayi [67], P. jishanensis [68] and cultivated P. rockii [52,69], but lower than that in wild P. rockii [70]. These findings suggest that genetic diversity varies significantly among different populations, with wild populations potentially harboring a greater degree of genetic variation. Overall, the SSRs and plant materials utilized in this study exhibited significant allelic diversity. Moreover, the disparities among different primers and plant materials could influence the results.

Additionally, the FIS for the population range from −0.95 to 0.82, with an average of 0.040, and 46 SSRs presented negative values. This result suggests a potential excess of heterozygotes in the population consisting of 123 P. ostii individuals, which may indicate a significant heterozygote advantage or outbreeding within the syngen. We propose that this phenomenon is closely associated with interspecific hybridization underlying horticultural tree peony evolution and self-incompatibility, which aligns with previous research findings [52,66]. Moreover, during the evolutionary process, the mutation sites were under the effect of positive selection pressure, and the influence of heterosis further exacerbated the complexity of the genome, thereby forming the highly heterozygous and redundant genome of P. ostii [27]. Therefore, our research not only confirmed that high heterozygous redundancy is a significant characteristic of P. ostii but also further supported the view that cultivated tree peonies have a complex hybrid origin [18,19].

4.2 LD and Genetic Structure

Evaluating the LD level in a population not only establishes a fundamental genetic basis for association studies but also guides the strategic selection of optimized analysis approaches. In the present study, the LD level between different SSR markers was relatively low. In general, outcrossing woody plants exhibit low LD levels because of their reproductive traits [71], and the characteristics of tree peony, as a typical cross-pollinated plant, are similar to those reported in other related studies [47,52,72]. We further speculate that extensive human intervention measures, such as introduction, selective breeding, and controlled pollination, may constitute some of the key factors contributing to the low LD level observed in P. ostii. However, for the association population in this study, the precise mechanism underlying the LD level remains unclear, primarily because the exact chromosomal distances of these loci have yet to be fully elucidated. Furthermore, low LD levels at a few SSR loci do not adequately represent the overall LD levels across the entire genome or intergenic regions [73]. Consequently, to comprehensively evaluate the LD level of a population, it is essential to employ multiple markers distributed throughout the genome for analysis, as different marker types may yield distinct insights because of their unique characteristics. In this study, we primarily employed a single-marker MLM for association analysis. This choice was motivated by (1) the observed low genome-wide LD in our population, which reduces the advantage of multilocus models that rely on extensive LD for haplotype construction; and (2) the goal of identifying individual, highly informative EST-SSR markers that could be directly applied in marker-assisted selection for their simplicity and cost-effectiveness. Beyond evaluating the LD level in the target population, a thorough examination of its genetic architecture is imperative. In actual studies, owing to the influence of various factors, it is nearly impossible to achieve a situation where there is no population structure at all. Existing studies have demonstrated that sample structure significantly influences genetic association analysis and has been recognized as a critical factor contributing to false associations [74]. Therefore, predicting population structure is among the primary conditions for conducting association analysis, which not only increases the accuracy of the analysis but also effectively prevents the occurrence of false positives. This research employed a triad of integrative approaches to characterize genetic architecture, revealing high consistency across methodological outcomes. Through STRUCTURE analysis, we subdivided the associated population into two primary subgroups, a finding that was strongly corroborated by both NJ tree and PCA analyses. This finding contrasts with a previous report suggesting the existence of eight subgroups in P. ostii [66]. We speculate that this discrepancy could be attributed to the different types of EST-SSR markers used and the varying sizes of the population. Despite accounting for subpopulation stratification within the target population, the general linear model (GLM) remains limited in fully mitigating spurious associations [75]. Conversely, MLM demonstrates superior efficacy in terms of suppressing errors. The joint incorporation of population structure (Q) and kinship (K) matrices markedly decreases error rates [76]. Additionally, applying FDR adjustment to all trait-associated p-values improves result reliability by effectively mitigating inflation.

4.3 Associations with Hundred-Grain Weight

SSR markers derived from candidate genes demonstrate greater genetic influence in modulating the transcriptional activity and functional roles of quantitative trait-linked genes [77]. Within the mapping population of P. ostii, 21 EST-SSRs exhibiting significant associations with hundred-grain weight were discovered in this study. Notably, SSRs from gene families such as AP2, YAB, and B3 exhibited pronounced genetic effects, suggesting that the target trait may be coregulated by multiple gene families and exhibit colocalization. These findings further corroborate the theory that quantitative traits are typically governed by the synergistic action of multiple minor-effect genes. Moreover, we found that certain SSR markers were significantly correlated with multiple traits. For example, P280 and PS27, two AP2 transcription factors, were significantly correlated not only with the hundred-grain weight but also with multiple fresh biomass traits in tree peony, including fruit weight, seed weight of a fruit, fruit yield per plant, and seed yield per plant [72]. This pattern might arise because of inherent significant correlations among these phenotypic traits and further highlights the characteristics of gene pleiotropy [78]. Consistent conclusions have also been reported in studies involving other woody plants [79]. These pleiotropic associations not only help identify key genomic regions but also highly important for achieving trait improvement through molecular marker-assisted breeding.

Association studies targeting complex quantitative traits in plants can identify numerous markers that are significantly correlated with these traits. However, these markers collectively explain only a limited proportion of the observed phenotypic variation [75]. In the current research, the mean phenotypic variance explained by EST-SSR markers for hundred-grain weight reached 18.61%, notably surpassing the explanatory rates documented in earlier works: 5.50% for floral traits [52] and 6.53% for fruit traits [47]. Furthermore, an in-depth analysis of the genetic regulatory mechanisms underlying quantitative traits in plants can provide a robust scientific foundation for conducting breeding work using significant association combinations. For instance, loci exhibiting overdominance effects identified through association analysis suggest that compared with homozygotes, heterozygotes may display superior phenotypic traits under certain conditions. In this research, PS55 and PS142 were significantly correlated with hundred-grain weight and presented an overdominance effect (with |d/a| values of 37.63 and 34.40, respectively). These findings indicated that individuals possessing the heterozygous alleles of PS55 and PS142 within the associated population might produce seeds with increased weight. Therefore, the integration of over-dominant genetic loci could strengthen heterosis performance. In summary, by integrating genotypes that display consistent effects, an effective strategy can be provided for the early screening of desired agricultural characteristics.

The repeated verification of genotype-phenotype association loci is highly important for association analysis, as it can effectively reduce the occurrence of false-positive associations. Among the 81 markers used in this study, 7 were consistent with those reported by Cui et al. [62]. Nevertheless, no identical association combinations were identified in the association analysis of hundred-grain weight between the two studies. This outcome likely stems from multiple causes, including limited sample size and multifaceted gene-environment interactions, diverse genetic backgrounds (e.g., differences between P. ostii and P. rockii), and gene-gene interactions. Consequently, some genuine associations might be challenging to replicate [79]. Further research should integrate diverse germplasm resource populations and combine phenotypic data collected across multiple years and various sites to perform comprehensive validation of EST-SSRs [80].

The fundamental strategy of genome-wide association studies relies on high-density genotyping of markers spanning the entire genome. Such an approach ensures that functional alleles are in linkage disequilibrium with at least one marker. The initial process step involves screening adequate polymorphic markers (e.g., SSRs/SNPs). For instance, genomic studies have demonstrated that ~130,000 markers achieve near-complete coverage in Arabidopsis thaliana (genome size: 125 Mb), which aligns with current sequencing standards [81]. Similarly, approximately 2 million markers are estimated to be required for effective coverage of the Vitis vinifera genome, which is approximately 475 Mb in size. Furthermore, characterizing the genetic diversity among different maize varieties may necessitate the use of as many as 10 to 15 million markers [82]. In the present study, we screened 81 SSR markers and applied them to association mapping analysis. Given the limited number of markers, the findings of this study do not constitute a comprehensive genome-wide association analysis. Nevertheless, this study represents the first association analysis conducted on the hundred-grain weight of P. ostii, providing a solid foundation for subsequent in-depth association analysis. As the cost of sequencing and genotyping technologies continues to decrease, more complex and refined association analyses with higher marker density will become increasingly feasible.

4.4 TFs Associated with Hundred-Grain Weight

P. ostii has emerged as a novel woody oil crop in China, and the hundred-grain weight serves as a critical factor influencing its yield. Studies have also demonstrated that TFs from the AP2/ERF gene family play a crucial role in increasing crop yield [83,84,85]. As a member of a plant-specific transcription factor superfamily, AP2/ERF critically regulates morphogenesis, stress adaptation, and metabolic homeostasis through conserved molecular mechanisms [86]. The present study revealed that TFs from AP2 not only strongly correlated with hundred-grain weight but also exhibited relatively high association frequency. In Arabidopsis and rice (Oryza sativa), AP2-domain TFs are conserved in terms of their ability to coordinate seed development and yield-related traits [83]. In this study, the AP2 family gene WRI1 (P280) was identified as a key regulator of P. ostii seed weight, as demonstrated by a significant correlation between the hundred-grain weight (|d/a| = 3.78) and overdominance genetic effects (|d/a| > 1.25). The observed dosage effects position P280 as a promising molecular marker for enhancing seed yield in this oil-producing peony species. In addition, WRI1 has been identified in both monocotyledonous and dicotyledonous plants, including P. ostii, Ricinus communis, Avena sativa, and Camelina sativa, clarifying its expression pattern during seed development and its key role in regulating lipid accumulation [53,87,88,89,90]. Related studies have demonstrated that WRI1 functions in regulating the expression of related genes in the glycolysis and fatty acid biosynthesis pathways, thereby significantly contributing to the biosynthesis of plant oils [91,92,93]. With oil contents ranging from 27–33%, the oil concentration and seed weight of tree peony seeds are directly correlated [31]. In addition to WRI1, another member of the AP2 gene family, ERF7 (PS27), was significantly correlated with the hundred-grain weight, demonstrating an explanatory power of 29.78%. Plant nitrogen metabolism critically regulates developmental processes and significantly influences crop yield. Previous research has indicated that ERF7 (PS27), GATA8 (PS118) and WRKY22 (PS131) participate in nitrogen metabolism by forming an interactive regulatory network. Among them, the synergistic effect of ERF7 (PS27) is reflected in its gene effect value (|d/a| = 0.26 < 0.50, indicating the presence of additive effects). Moreover, previous studies have demonstrated that synergistic interactions among related genes can effectively delay plant senescence and substantially increase crop yield [94,95].

In this study, NGAL3 (PS98) was also significantly associated with hundred-grain weight. Members of the RAV subfamily within the B3 family of transcription factors, specifically NGALs in Arabidopsis, play crucial roles in regulating plant growth and organ development [96,97]. Among these genes, NGAL2 and NGAL3 cooperate with KLU in a shared regulatory pathway that controls seed size. Specifically, NGAL2 binds to the CACTTG motif within KLU promoter, suppressing its expression. This downregulation consequently impedes maternal outer integument cell division and reduces cellular proliferation in these structures, ultimately yielding relatively small seeds [57]. Additionally, in addition to regulating seed and flower development, NGAL3 can also inhibit the expression of CUC2 in a miR164-independent manner, consequently influencing leaf margin development [98,99,100].

In addition to the above three TFs, YAB1 (PS95–96) was significantly associated with the hundred-grain weight of P. ostii, explaining 33.04% of the variation. The Arabidopsis thaliana genome encodes six YABs that are phylogenetically clustered into five distinct subfamilies. Among them, AtFIL/AtYAB1, AtYAB2, AtYAB3 and AtYAB5 jointly regulate the development of organs such as flowers and leaves, whereas AtCRC and AtINO mainly participate in the formation of floral organs, such as carpels, nectaries and ovules [101]. In Setaria italica, the overexpression of SiDL (one of the SiYABs) further confirmed these functions in Arabidopsis, leading to leaf curling, delayed flowering, reduced seed size, and negative regulation of the response to salt stress [102]. Rice (Oryza sativa) possesses eight YABs, of which OsYAB1 executes feedback control in GA biosynthesis through the transcriptional regulation of two pivotal enzymes: GA3ox2 and GA20ox2. Importantly, ectopic expression assays revealed the dual functionality of this gene—enhancing carpel and stamen proliferation—while confirming its regulatory role in floral morphogenesis [103]. Grapes (Vitis vinifera) have seven YAB genes, with VvYABBY4 experimentally validated as a key regulator during seed morphogenesis. Collectively, these genes regulate the development and maturation of seeds, thereby influencing the reproductive traits and yield of grapes [104,105]. This studies have shown that the YAB gene family serves as a critical regulator in the flowering process of plants, regulating the normal development of organs such as leaves and flowers, as well as the vitality of meristems [106]. In the present study, YAB1 was significantly associated with the hundred-grain weight of P. ostii (with an explanatory rate of 33.04%). Therefore, we speculate that YAB1 ultimately affects the hundred-grain weight of tree peony seeds by regulating their development and maturation.

In summary, the gene loci of WRI1 (P280), ERF7 (PS27), NGAL3 (PS98) and YAB1 (PS95–96) identified in this study may be key candidate genes that regulate the weight of tree peony seeds; these genes have important research value and warrant further in-depth exploration. Future perspectives: To functionally validate the candidate genes identified in this study, a clear research roadmap is warranted. The immediate next steps include (1) performing qRT-PCR analysis to examine the spatiotemporal expression patterns of these genes throughout seed development in P. ostii accessions with different seed weights; (2) performing allele-specific expression analysis for polymorphic sites within these genes; and (3) implementing genetic transformation (e.g., overexpression or CRISPR-Cas9 knockout) in model plants such as Arabidopsis or tobacco to directly assess their effects on seed size and weight.

5 Conclusions

In conclusion, this study successfully delineated the genetic architecture underlying the hundred-grain weight in the emerging woody oil crop P. ostii. By constructing a genome-wide association map and employing polymorphic EST-SSR markers, we identified two distinct genetic subgroups within the population and pinpointed 21 marker–trait associations, which prominently involved several key transcription factor families (including AP2, TCP, and YABBY). Further analysis revealed WRI1, ERF7, NGAL3, and YAB1 as the most promising candidate genes for regulating seed weight. These findings extend beyond mere association by offering concrete genetic targets and a robust theoretical framework for the molecular breeding of P. ostii. While these findings provide a critical foundation, future functional validation of these candidate genes through techniques such as gene editing or transgenics will be essential to confirm their roles and ultimately translate this knowledge into tangible yield improvements in woody oil crops.