The recombinant inbred line (RIL_F₇) population was derived from two parental lines, Y3 and K326 through the single-seed descent method. This population contains 271 genotypes and was planted at Shilin in a completely randomized design in 2020, 2021, and 2022 years, which were treated as three distinct environments. The evaluated traits for each genotype included total nitrogen % (TN), potassium % (POT), chlorine % (CHL), reducing sugar % (RS), total plant alkali % (TPA), total sugar % (TS), starch % (STA), chlorogenic acid mg/g (CHA), rutin mg/g (RU), fructose % (FRUC), xanthophyll μg/g (XAN), beta-carotene μg/g (BCA), citric acid mg/g (CA), petroleum ether % (PE), cellulose % (CE), the difference between two sugars % (DS), total amino acids (TAA), aspartic acid (APA), phenylalanine (PA), glutamine (GL), and protein % (PRO). These traits were quantified using high-performance liquid chromatography (HPLC) following the procedures detailed in (Julio et al., 2006; Pang et al., 2006; Jing et al., 2016; Wang P. et al., 2024). Eighteen chemical traits were evaluated in E1 (2020-SL), nine traits in E2 (2021-SL), and eleven traits in E3 (2022-SL).

We estimated variance components using the following mixed linear model (Tong et al., 2025).

In this model, Y k h represents the phenotypic value of the k-th genotype in the h-th environment; μ represents the population mean; g k indicates the genotypic value of the k-th genotype, g k ~ N ( 0 ,   σ g 2 ) ; e h ~ N ( 0 , σ e 2 ) , denotes the effect of the h-th environment; ϵ k h ~ N ( 0 , σ ϵ 2 ) , denotes the residual effect of the k-th genotype in the h-th environment. To estimate variance components ( σ ^ g 2 , σ ^ e 2 , σ ^ ϵ 2 ), mmer module of the Sommer R package was utilized and genotypic values were predicted using the best linear unbiased prediction (BLUP) method. Broad sense heritability was calculated using the formula H 2 = σ ^ g 2 / ( σ ^ g 2 + σ ^ ϵ 2 ) , here σ ^ g 2 represents the genotypic variance and σ ^ ϵ 2 represents residual variances (Tong et al., 2025). Additionally, the pearson correlation coefficient between traits was computed using the rcorr function of the Hmisc R package. Phenotypic correlation coefficients were calculated separately for each environment, while the genetic correlation coefficients were derived from predicted genotypic values using the BLUP method. Finally, the relationship between traits was visualized using the corrplot package in R.

A total of 274 samples, including two parental lines, one F₁ hybrid, and 271 F₇ individuals, were genotyped using the BIGSEQ-500 platform, following the protocol outlined in our previous study (Tong et al., 2023). High-quality reads were filtered and aligned to the reference genome Nitab4.5 (Edwards et al., 2017) using bioinformatics tools. SNPs and InDels were called out using GATK, with stringent quality control filters. From these data, 46,324 bin markers were constructed and used to develop a high-density linkage map. This map spans a total genetic distance of 3334.88 cM across 24 linkage groups (LGs), with an average marker interval of 0.469 cM ( Supplementary Table S1 ).

We employed a full QTL model to investigate the genetic architecture of complex traits across multi-environment field experiments. This model incorporates the individual additive genetic effect (a) of each QTL, the additive-by-additive epistatic effect of each QTL pair (aa), and their interaction with the environments (ae and aae). We assume ‘s’ is the number of segregating QTLs and ‘t’ denotes the number of QTL pairs exhibiting epistasis. Then, the phenotypic value of the k-th genotype in the h-th environment can be described by the following mixed linear model (Tong et al., 2023):

Where, μ is the population mean; a_i is the additive effect of the i-th QTL with coefficient x i k which is treated as a fixed effect and takes values 1 and − 1 for QQ and qq genotypes of QTL, respectively. Similarly, a a i j is the additive-by-additive epistatic effect of the i-th and the j-th QTL with coefficient x i k x j k as a fixed effect; e h is the random effect of the h-th environment, e h ~   ( 0 , σ e 2 ) ; a e h i indicates the additive by environment interaction effect of the i-th QTL and the h-th environment with coefficient u h i k ( = x i k ) , a e h i ~ ( 0 , σ a e i 2 ) ; a a e h i j refer to the interaction effect of the h-th environment with a a i j , with coefficient u h i j k ( = x i k x j k ) , a a e h i j ~ ( 0 ,   σ a a e i j 2 ) ; ϵ k h is the random residual effect, ϵ k h ~   ( 0 , σ ϵ 2 ) .

QTL analysis was conducted using QTLNetwork 2.0 (Yang et al., 2008), specifically designed for mixed-linear-model-based composite interval mapping (MCIM). Both one-dimensional (1D) and two-dimensional (2D) genome-wide scans were conducted at a walking speed of 1 cM, To control experiment wise type 1 error rate, a critical F-value based on Henderson method III which was determined by permutation testing 1000 times for each tested locus at a significance level of 0.05. The full QTL model was employed to estimate and test the QTL effects and their significance using the Markov Chain Monte Carlo (MCMC) method. Finally, the distribution of QTLs across linkage groups ( Figure 1 , Figure 2 ) was visualized using the LinkageMapView package in R.

Each QTL region was defined by the two flanking bin markers in the genetic linkage map. The sequences of these markers were aligned to the Nitab4.5 reference genome using Burrow Wheeler Aligner with the mem algorithm (Li and Durbin, 2009) and the genes were extracted using intersect function in BEDTools (Quinlan and Hall, 2010). Variants including SNPs and Indels were annotated with SnpEff (Cingolani et al., 2012) and those predicted to have moderate to high impact on protein function were retained for further analysis. These variants were subsequently validated through single marker association analysis using PLINK (Purcell et al., 2007), applying a significance threshold of p < 0.05. For functional enrichment analysis, the protein sequence of Nitab4.5 reference genome was uploaded to eggnog-mapper website. Finally, Gene Ontology (GO) and KEGG pathway enrichment analysis were performed using clusterprofiler R package (Yu et al., 2012). Putative candidate genes were functionally characterized using the BLASTp module of NCBI against the non-redundant (nr) protein database.

For multi-environment phenotypes, variance components and broad-sense heritability were estimated ( Table 1 ). The heritability of traits TS, RS, TN, and POT was 0, indicating that these traits were predominantly influenced by environmental and error variance. In contrast, traits DS and TPA were significantly influenced by genetic factors, as reflected by their heritability of 75.71% and 44.83%, respectively. The majority of traits exhibited statistically significant but negative correlation values (α=0.05) ( Figure 1 ). Significant phenotypic correlation was observed between TS and RS across all environments, followed by APA and PA in E1 ( Figure 1A ). The correlation between PA and GL in E1 exhibited a similar trend, as did TN and PRO in E2 and E3, respectively ( Figures 1B, C ). In contrast, a consistently large negative correlation was observed between TS and TN across all environments, though the magnitude varied. The phenotypic correlation between CHL and STA was positive and significant (α=0.05) in E1 but negative and non-significant in E2 and E3. Overall, the correlation pattern varied slightly across environments but remained consistent with general trends. The highest genetic correlation coefficient of 0.33 was observed between POT and PRO ( Figure 1D ), accompanied by a substantial phenotypic correlation in E2 and E3 ( Figures 1B, C ). Additionally, the traits TS, RS, and TN exhibited no genetic correlation, demonstrating that environmental influences drove phenotypic variation rather than shared genetic architecture.

We identified 18 QTLs associated with 12 traits that exhibited significant individual additive effects across 5 linkage groups ( Figure 2 ) and three pairs of epistatic QTLs distributed over 6 linkage groups for two traits ( Figure 3 ). Two epistatic QTL pairs were detected for chlorogenic acid (CHA) and one QTL pair for rutin (RU). LG15 contained the highest number of QTLs with individual additive effects (10 QTLs), followed by LG06 with 5 QTLs. LG11, LG16, and LG24 each carried one QTL.

A total of 18 QTLs with additive (a) effects were detected for 12 traits ( Table 2 ). Out of these 18 QTLs, 11 QTLs contributed positive additive effects, while 7 QTLs exhibited negative additive effects, indicating the complex genetic architecture ( Supplementary Figure S1 ). The heritability of each QTL explained a percentage of the phenotypic variation, ranging from 2.21% to 20.05%. The majority of the QTLs had small additive effects and lower heritability, thus regarded as minor-effect QTLs. The QTLs qAPA16-247 ( h a 2 = 20.05 % ) , qGL15-18 ( h a 2 = 18.85 % ) , qAPA15-250 ( h a 2 = 15.75 % ) , qPA15-249 ( h a 2 = 12.2 % ) , and qTAA15-248 ( h a 2 = 11.16 % ) were identified as major-effect QTLs due to their significant and larger contribution to phenotypic variation. However, qTAA15-248, qTAA15-169, qGL15-18, qAPA15-250, qPA15-249, and qPA15-18 exhibited substantial additive effects on their respective traits and were located on LG15. For breeding purposes, one of the main goal is to identify QTLs that express stably across environments with minimal or non-significant QTL-environment interactions. In our study, all QTLs exhibited no interaction effects, indicating they can be utilized in breeding new variety for most general environment. Notably, qPA15-18 and qGL15-18 were located in a same marker interval ranged by SNP_0010516_37142 and SNP_0198499_1362. The co-localizations indicate there may be gene which take pleiotropic effect on PA and GL. This hypothesis was further supported by the substantial and reasonably high estimated correlation coefficient of 0.87 between PA and GL ( Figure 1A ).

A two-dimensional (2D) genome wide scan detected three epistatic QTL pairs associated with chlorogenic acid and rutin across LG1/LG16, LG2/LG11, and LG5/LG7 ( Table 3 ). All these QTL pairs exhibited minimal additive-additive epistatic effects. Furthermore, these epistatic QTLs exhibited no significant additive-additive epistasis by environment interaction effects. Each QTL pair explained less than 3% of the overall phenotypic variation. Notably, the QTL pair qRU5-182/qRU7-131 accounted for greater heritability ( h a a 2 = 2.37 % ) than qCHA1-216/qCHA16-21 ( h a a 2 = 1.45 % ) and qCHA2-281/qCHA11-56 ( h a a 2 = 1.29 % ) . Interestingly, the interaction effects of qCHA2-281/qCHA11-56 and qRU5-182/qRU7-131 were -0.26 and -0.27, respectively ( Supplementary Figure S2 ). This indicated that the genotype of two QTL from same parent will reduce the trait value, in contrast, the genotype from different parent will increase the trait values.

Through comparative mapping with the Nitab4.5 reference genome, the QTLs were mapped onto ten chromosomes, namely Nt05, Nt06, Nt08, Nt10, Nt12, Nt16, Nt18, Nt21, Nt22 and Nt24. A total of 477 genes were identified from these QTL regions ( Supplementary Table S2 ). 56,983 variants within these genic regions were identified and annotated through the SnpEff tool. Wherein, 395 variants in 99 genes were determined to have moderate to high impacts on the protein level and were filtered for further analysis ( Supplementary Table S3 ). These 99 genes were subjected to genome-wide association study using single marker association analysis. This analysis revealed 205 variants in 66 genes demonstrating significant association with multiple phenotypes across all the environments ( Figure 4 ), ( Supplementary Table S4 ). Among the 66 genes, 20 were significantly enriched in six GO biological processes, eight in GO cellular components, and one in GO molecular functions, while three genes were enriched in four KEGG pathways ( Figure 5 ). Based on these results, Nt08g00266, Nt16g00236, and Nt22g03479 were predicted as candidate genes. GWAS analysis further revealed that Nt08g00266 was significantly associated with total sugar TS and STA, Nt16g00236 with TPA, and Nt22g03479 with TS and RS. The functions of these candidate genes were retrieved from BLASTp. This analysis revealed that Nt16g00236 encodes a mitogen-activated protein kinase (MAPK), Nt22g03479 encode scopoletin glucosyltransferase, and Nt08g00266 encodes a MYC2-like transcription factor.