An interspecific F₁ hybrid population derived from a cross of Daye (L. chinense) as the female parent and ZKLC1 (L. barbarum) as the male parent (Rehman et al., 2020) was used in the present study. Of the 450 hybrid progenies, 305 were chosen at random in September 2018 and were used to create a genetic linkage map and QTL mapping study of fruit size and quality-related characteristics. The botanical features of the female parent Daye, also referred to as Chinese goji (L. chinense), include ovate to lanceolate leaves, shallow to dark purple flowers, an oblong fruit shape, fruit color ranging from dark red to orange, fruit weight (0.38 g), fruit length (1.08 cm), fruit diameter (0.75 cm), and seeds that are generally broader, rounded, and yellowish to yellow-brown in color. The male parent, ZKLC1 (L. barbarum), has round fruit with a diameter of 1.28 cm, weight of 1.05 g, length of 1.79 cm, and color of light-yellow seeds that are semi-spherical to flattish. It also has lanceolate to oblong leaves and vivid to royal purple flowers (Rehman et al., 2020). Fully mature fruits of this interspecific hybrid population in the summer season were used as the experimental materials. Fruits of 305 F₁ individuals and their parents were collected (50 fruits from each tree of similar size and maturity were selected) and stored in a low-temperature refrigerator at 4°C.

Goji fruit tasting evaluation was performed in the laboratory at room temperature at the Northwest China Bio-agricultural Center (38°28′05″ N, 106°16′23″ E), Yinchuan City, Ningxia Hui Autonomous Region, China. Thirty bitter and sweet sensory evaluation forms were prepared, including the name of the participant, gender, age, and a scale for bitterness from 1 (Tasteless) to 7 (extremely bitter) and similarly for sweetness from 1 (Tasteless) to 7 (extremely sweet). In addition, we arranged materials such as mineral water, tasteless cookies, and spittoons for each evaluator. For the taste survey, we invited 30 people (male and female) aged 20 to 60 to evaluate the goji fruit taste from 7:30 to 11:30 AM and from 14:00 to 18:00 PM and completed it on the same day. The taste assessment process was as follows. First, a staff member selected fruit samples from an F₁ individual and distributed one of its fruits to each individual. The taster then chews the fruit completely in his or her mouth to fill it, spits it out into a spittoon, and rates the sweetness and bitterness within one minute. Finally, gargle with mineral water and get ready to taste the fruit of the next individual and so on. Each person from the panel of 30 people had to rate the goji fruit taste (sweetness and bitterness) for all 305 samples on the same day. Consecutively, red goji fruit taste data were evaluated for three years (2020–2022) during the summer season by the same panel of 30 people each year.

Five fruits of the same maturity were chosen for each F₁ individual and their parents, and the soluble solid content (SSC) (°brix percentage) was measured using a digital pocket refractometer ATAGO PAL-1 (Atago Co. Ltd., Tokyo, Japan). The SSC data were collected consecutively for four years (2019–2022). The results are expressed as °brix percentage. Whereas fruit firmness (FF) was also determined for the five mature fruit repeats per hybrid progeny using the fruit hardness tester equipment GY-4 (Zhejiang Top Instrument Co., Ltd. Hangzhou, China) (Rehman et al., 2020), and data were collected consecutively for three years (2019–2021). The operation of the instrument involves positioning each fruit horizontally on the base and applying gentle, uniform pressure with the pressure head of the penetrometer into the flesh of the fruit. The sensor then displays the reading on the digital panel, and the average value of five consecutive measurements is obtained for each hybrid progeny. Following completion of each set of readings, the scale was reset to zero. Here, we used one-year SSC and FF (2019) data from our previous study (Rehman et al., 2020), and newly collected SSC (2020–2022) and FF (2020–2021) from the F₁ interspecific hybrid population of red goji berries to re-analyze and update the QTL mapping results based on the reference genome of goji berries (Cao et al., 2021).

The morphological data of 305 F₁ individuals of fruit size-related traits, such as fruit weight (FW) in grams, fruit length (FL), fruit diameter (FD) in centimeters, and fruit shape index (FSI) (determined by dividing fruit length by diameter), were consecutively collected for three years (2018–2020) and statistically analyzed. Data were collected as previously described (Rehman et al., 2020). In the current study, we used two years of data (2018–2019) of four traits (FW, FL, FD, and FSI) from our previous study (Rehman et al., 2020), along with another year (2020) of data to re-analyze and update the QTL mapping results using the goji reference genome (Cao et al., 2021).

Analysis of variance (ANOVA) for three or four individual years, Pearson correlation, and general descriptive statistical analysis, including frequency distribution, were carried out using Origin Pro v. 2024, OriginLab Corporation, Northampton, MA, USA. (https://www.originlab.com/).

The 305 hybrid progenies and parents’ young leaves were collected and dried using a desiccant (silica gel) before DNA extraction. Using a slight modification of the manufacturer’s instructions, DNA from each leaf was extracted using the Plant Genomic DNA Kit (TIANGEN BIOTECH (BEIJING), CHINA, CO., LTD.). The quantity and quality of extracted DNA samples were assessed using a Thermo Fisher Scientific NanoDrop ND-2000 spectrophotometer and 1% agarose gel electrophoresis, respectively (Rehman et al., 2020). Resequencing was performed according to a standard protocol established by Illumina. The genomic DNA of each sample was subjected to library construction once it passed the detection. Those that passed quality control were subjected to the Illumina HiSeq platform for further sequencing. To ensure data quality, raw reads were filtered to obtain paired-end sequencing data following the steps of data filtering, such as the removal of joint sequences contained in reads, using fastp software v. 0.23.0 (Chen et al., 2018); the low-quality bases among reads were removed (the average mass number was calculated using a sliding window of 4bp, and all bases were removed below 15); reads must be greater than 50bp in length. MEM algorithm of the BWA software v. 0.7.15-r1140 (Li, 2013), was first used to compare and align the filtered reads with in-house data along with the L. barbarum reference genome with a genome size of 1.68 Gb (Cao et al., 2021) (only the paired-end reads that matched both ends were considered suitable for genome comparison). The SAM tools v. 1.3.1 (Li et al., 2009) was used to exclude duplicate reads from the aligned reads. Genomic variations in each accession were determined using the Genome Analysis Toolkit (GATK) v. 3.7 software’s Haplotype Caller module and GVCF model (McKenna et al., 2010). All GVCF files were merged. Using the Haplotype Caller module, we filtered variants including SNPs and InDels based on strict filtering criteria. The selection criteria were as follows: (a) at least one parent was heterozygous, (b) depth for parents ≥5, (c) depth for individual ≥5, (d) the proportion of missing data in the hybrid progeny did not exceed 25%, and (e) a minor allele frequency (MAF) of less than 10% (Catchen et al., 2011). The SNPs and InDels that were identified were combined and turned into bin markers (Qi et al., 2021). These bin markers were then divided into separate genomic regions and genes, and further annotated using ANNOVAR v. 2016 Feb 1st (Wang et al., 2021) to predict gene function. The region of the mutation site in the genome (intergenic region, gene region, CDS region, etc.) and the impact of the mutation (synonymous and non-synonymous mutations, etc.) were determined based on the location of the mutation site in the reference genome and gene location information on the reference genome.

The construction of the linkage map harboring filtered bin markers was facilitated by the use of Lep-MAP3 software with default parameters, which is based on the maximum likelihood method as described by Rastas (2017). To ensure the precision of the genetic map, bins less than 5kb were initially filtered out, and bins that deviated from the expected Mendelian segregation with a significance level of P < 0.001 were excluded from the analysis. The total length of the genetic maps was estimated using the Kosambi mapping function initially proposed by Kosambi (1943). This function is commonly employed to convert recombination frequencies to map distances in centiMorgans (cM). The Order-Marker2 module was used to sort the markers within each linkage group, calculate the genetic distance of each LG, and verify the continuity of the generated maps. A suite of Python MadMapper scripts was used to evaluate the quality of the labeled bin markers on the genetic linkage maps (https://cgpdb.ucdavis.edu/XLinkage/MadMapper/). To assess the quality of the genetic linkage map, a heat map of recombination rates between markers, haplotype map for each hybrid progeny, and collinearity map of the genetic map and genome were generated using CIRCOS v. 0.66 (Krzywinski et al., 2009). GACD software v. 1.2.13 (Zhang et al., 2015), was used for QTL or genomic region identification based on the inclusive composite interval mapping (ICIM) model (Broman and Speed, 2002; Li et al., 2008). The significant threshold was manually set at 2.5 for LOD scores, with 1.0 cM as the walking speed for all QTLs and 0.001 as the PIN value (Wang, 2009; Zhou et al., 2024). QTLs of the eight traits associated with bitter taste (BT), sweet taste (ST), soluble solid content (°brix percentage) (SSC), fruit firmness (FF), fruit weight (FW), fruit length (FL), fruit diameter (FD), and fruit shape index (FSI) were detected across three or four environments over a period of three- or four-years data. A major QTL was considered a QTL with an LOD value > 3 and a phenotypic variance of > 10% across several environments, whereas a stable QTL was found in at least two or three separate environments (Fan et al., 2017; Che et al., 2018). If the confidence intervals overlapped, the QTLs for identical traits found in the other environments were deemed to be the same. The possible locations of the QTLs were described based on the LOD peak sites and surrounding regions. According to the computed results by GACD, which have been thoroughly explained by Nzuki et al. (2017), the additive (a) and dominance (d) effects were calculated using the Muchero formulation (Muchero et al., 2013). The QTL mode of action was determined by dividing dominance by the absolute additive value (d/|a|) ratio. A ratio greater than one was considered over-dominant, a ratio between 0 and 1 as partial dominance, and a ratio less than 1 as under-dominance (Hua et al., 2003).

Fruit quality and fruit size/weight related traits, including ST, BT, SSC, FF, FW, FL, FD, and FSI, were evaluated from 305 F₁ hybrid progeny of interspecific population based on three or four individual years. Analysis of variance (ANOVA) of the studied traits showed significant differences (P < 0.05, P < 0.01) among F₁ progenies in different years ( Supplementary Table 1 ). Descriptive statistical analysis showed that the coefficient of variation (CV%) ranged from 16 to 22 for ST, 17 to 21 for BT in three years (2020–2022), 8–11 for SSC in four years (2019–2022), 19–22 for FF in three years (2019–2021), 21–22 for FW, 9–12 for FL, 7–10 for FD, 8–9 for FSI in three years (2018–2020) ( Supplementary Table 2 ). The frequency distribution histogram and box chart of all eight traits for individual years showed a normal distribution and significant differences among the different years, respectively ( Figures 2A–H ; Supplementary Figures 1A–H ). Normality tests were performed using the Kolmogorov-Smirnov (K-S) goodness of fit based on the absolute distance between the cumulative distribution and values ranged from 0.02 to 1, suggesting a positive normal distribution of all traits evaluated under different years except BT21, BT22, ST22, SSC22, and FW19 ( Supplementary Table 3 ). Correlation analysis (P < 0.05 or P < 0.01) revealed an extremely significant positive association detected in a few comparisons among BT and FW, FD; ST and SSC, FW, FD; SSC and FD; FF and FW, FL, FD, FSI; and fruit size related traits such as fruit weight, fruit length, fruit diameter and fruit shape index showed highly positive significant correlation among each other. A highly significant negative correlation was also detected between BT and ST, FW, FD; SSC and FSI, FW, FL, FD; and FD and FSI ( Supplementary Table 4 ; Figure 3 ).

Resequencing was performed for raw read data of the two parents and 305 F₁ individuals using in-house data and the L. barbarum reference genome (Cao et al., 2021). As shown in Supplementary Table 5 , approximately 36.27 and 29.67 million cleaned reads from parents (DY and ZKLC1) and an average of 9.46 million cleaned reads from F₁ individuals including 2.9 billion total clean reads, were obtained. The average GC content was 39.85% with Q20 and Q30 scores ranged from 98.91 to 96.54%. The average sequencing coverage and depth of DY, ZKLC1, and F₁ individuals were 5.31%, 49.33; 8.09%, 43.93; and 4.61%, 17.57; respectively. Of these reads, 50.62% from the female parent (DY), 82.19% from the male parent (ZKLC1), and 63.96% from F₁ individual plants were uniquely mapped onto the reference genome of L. barbarum and employed for SNP calling ( Supplementary Table 5 ). Finally, 3,058 bins containing total 35,331 polymorphic markers between the two parents were obtained, comprising 33,212 total SNPs with SNP density of 239.01 SNP/Mb, and 2,119 InDels with InDel density of 15.38 InDel/Mb. SNPs and InDels were distributed across the genome, as illustrated in Supplementary ( Supplementary Figure 2 ). Detailed information on these SNPs is provided in Supplementary ( Supplementary Table 6 ). Most annotated SNPs (79.5%) were detected in intergenic areas, and 3.9% were located in the exonic region. Of these, 54.3% were nonsynonymous SNPs. Compared to SNPs, intergenic regions accounted for over half (65.1%) of the annotations for InDels, with exonic regions accounting for 1.6% of the total. Frameshift mutations resulted in 67.6% exon InDels ( Figure 4 ).

Goji bin linkage map comprising 12 chromosomes were constructed harboring 3,058 bins containing 35,331 polymorphic markers, that varied from 197 to 377 bins and spanned a coverage ranged from 56.22 to 109.33 cM. The integrated genetic map spans a total length of 966.42 cM comprising 847.73 cM and 1065.89 cM lengths of female parent and male parent, respectively. The mean interval between bins were 0.03 cM with an average bin interval of 0.319 cM, and an average of maximum distance between the bins were 4.223 cM ( Supplementary Figure 3 ). The longest chromosome was Ch06 (Lba06) spanning a length of 109.33 cM and comprising 2,947 markers in 377 bins. Ch02 (Lba02), on the other hand, was the smallest chromosome, spanning 56.22 cM and having 197 bins with 5,823 markers ( Supplementary Table 7 ). The 1,121 bins on the female genetic map had an average maximum bin interval of 7.363 cM and an average spacing between the markers and bins of 0.027 and 0.786 cM, respectively. In contrast, the male genetic map had 1,966 bins, with an average maximum bin interval of 2.616 cM and an average spacing between the markers and bins of 0.033 and 0.551 cM, respectively ( Supplementary Table 7 ). The quality of the constructed consensus genetic map was evaluated using haplotype maps ( Supplementary Figure 4 ) and heatmaps ( Supplementary Figure 5 ) of 305 F₁ individuals based on 12 linkage groups. The results indicated that the bin genetic map was constructed with high quality and good collinearity corresponding to goji genome. To assess collinearity between the genetic map and the genome, 3,058 bins were mapped to the reference genome of the goji berries. For all 12 LGs, correlations between the genetic and physical maps were nearly linear ( Supplementary Figure 6 ). The order of the majority of bins in the linkage map matched that of the corresponding chromosomes in the physical map of the goji berry genome. Supplementary Table 8 presents the correlation coefficients for each LG. For all LGs, the average correlation coefficient between physical and genetic positions was 0.86. These LGs exhibited a high degree of genetic collinearity with the physical map, as demonstrated in our findings. Significant collinearity indicated that the markers adequately covered the genome of goji berries and appropriately covered the 12 chromosomes.

Using the ultra-dense genetic map of goji berries, a number of QTLs were plotted for physio-chemical and morphological traits based on three or four consecutive years of data. QTL mapping analysis was performed for BT, ST SSC, FF, FW, FL, FD, and FSI using the ICIM model in GACD v. 1.2.13. Most QTLs associated with the studied traits were located on distinct chromosomes, including chromosomes 3, 4, 5, 6, 7, 8, 9, and 12 ( Figure 5 ; Supplementary Table 9 ). A total of Ninety-nine QTLs targeting eight traits were found on separate chromosomes in the individual years ( Supplementary Table 9 ). Ten major effect and colocalized QTLs were detected on chromosomes 1, 3, 5, 6, 7, and 8, whereas stable QTLs targeting different traits were detected across two different environments ( Figure 5 ; Table 1 ). QTLs detected corresponding to each trait across three or four environments are described as follows:

Bitter Taste (BT): Seventeen QTLs targeting BT were detected on different chromosomes ( Supplementary Table 9 ), including two major QTLs (qBT-7.1, qBT-7.2) explaining the total phenotypic variance of 16.95%, 16.13%, and LOD scores of 29.25, 26.87, and positioned on chromosome 7 on the map with 20.5–21.5, 29.5–31.5, and flanking the marker interval (Marker13033977-Marker13326504; Marker19870129-Marker19934823) ( Figure 5 ; Table 1 ). Five loci (qBT-5.1, qBT-8.2, qBT-10.1, qBT-10.2, and qBT-6.1) targeting BT overlapped with (qFL-5.4-qFL-5.5, qBT-8.3, qST-10.3, qFL-10.4, and qFW-6.8), harboring a maximum PVE of 5.22% with 7.75 LOD and coincided confidence intervals ( Figure 5 ; Supplementary Table 9 ). The major QTLs showed positive and negative additive effects, indicating a large contribution from parents with partial dominance effects of the QTLs ( Table 1 ).

Sweet Taste (ST): Twelve ST QTLs were detected on different chromosomes ( Supplementary Table 9 ), including one major-effect (qST-12.1) and five minor-effect (qST-6.2, qST-10.3, qST-11.1, qST-11.2, and qST-12.2) QTLs were detected with the highest LOD score of 14.12 and PVE of 10.34% ( Figure 5 ; Table 1 ; Supplementary Table 9 ). The additive effect of major loci showed a positive contribution from both parents, with the partial-dominance effect of the QTL ( Table 1 ).

Soluble Solid Content (SSC) (°brix percentage): Sixteen genomic regions corresponding to SSC were located on different chromosomes ( Supplementary Table 9 ), including qSSC-6.4, a major-effect locus, with an LOD score of 15.24 and 11.50% variance. Additionally, 11 QTLs (qSSC-4.3, qFL-4.5; qSSC-6.5, qFL-6.10; qSSC-12.3, qSSC-12.4; qSSC-7.5, qFF-7.10; qSSC-7.6, qFF-7.8, qFF-7.11) were detected as overlapping loci corresponding to different traits with a maximum LOD of 7.25 and 6.22% phenotypic variance ( Figure 5 ; Supplementary Table 9 ).

Fruit Firmness (FF): A total of 10 FF QTLs were detected on different chromosomes with LOD ranged from 2.52–7.53 and PVE was 1.27–7.08% across the three environments ( Supplementary Table 9 ). A stable locus (qFF-7.8-qFF-7.11) was identified with an LOD score of 5.90 and phenotypic variance of 3.38 across two environments and showed negative additive effect of parents with an over-dominance effect of QTL ( Figure 5 , Table 1 ). qFF-12.5 was found to overlap with qFW-12.7 harboring maximum LOD of 7.79 and phenotypic variance 5% ( Supplementary Table 9 ).

Fruit size related traits: Thirteen FW QTLs were identified ( Supplementary Table 9 ), including three major-effect QTLs (qFW-1.5, qFD-1.8; qFW-3.6, qFD-3.7; qFW-6.6, qFW-6.9, qFL-6.13) accounting for LOD scores of 11.42, 9.01, and 10.35, and PVE of 12.85%, 11.27%, and 12.23%, respectively, and flanking the marker intervals (M3398563-M4171036; Marker127088931-Marker127359380; M127994636-M127996736), respectively ( Figure 5 ). These three major loci exhibited positive and negative effects on parents, indicating partial and over-dominance effects of QTLs, respectively ( Table 1 ). A stable locus (qFW-6.6, qFW-6.9, qFL-6.13) also overlapped with one colocalized locus (qFL-6.14, qFD-6.15) across the two environments, along with similar flanking marker and confidence intervals ( Figure 5 ; Table 1 ). The additive effect of the stable QTL was negative for parents, indicating over-dominance effects of the QTL ( Table 1 ). For FL, 16 loci were detected on different chromosomes ( Supplementary Table 9 ), including two major-effect and co-localized QTLs (qFL-8.7, qFD-8.9; qFL-8.8, qFD-8.10), accounting for LOD scores of 12.70, 26.63, and PVE of 11.16% and 23.05% and flanking the marker intervals (Marker30212799-Marker30212822; Marker33160675-Marker33394294), exhibited negative and positive additive effects with partial dominance effects of QTL ( Figure 5 , Table 1 ). qFL-2.3, and qFL-5.4, were observed to overlap with qFSI-2.4, and qFL-5.5, with a maximum LOD score of 4.59 and a variance of 2.34 ( Figure 5 ; Supplementary Table 9 ). Eight FD QTLs were detected on various chromosomes ( Supplementary Table 9 ), including qFD-8.8, a major-effect locus with an LOD of 27.08, and phenotypic variance of 14.25% ( Figure 5 , Table 1 ). Seven FSI QTLs were detected on different chromosomes ( Supplementary Table 9 ), including one major QTL (qFSI-5.7-qFSI-5.10) and two stable QTLs (qFSI-5.6-qFSI-5.9; qFSI-5.7-qFSI-5.10) with LOD values of 8.41 and 18.99, and phenotypic variance of 6.07% and 14.06%, respectively, flanking the marker interval (Marker2562548-Marker3019144; Marker6172058-Marker6899663) ( Figure 5 , Table 1 ). The additive effects of major and stable QTLs were negative and positive for parents, respectively, and showed partial and under-dominance effects of QTLs ( Table 1 ). The stable and overlapping genomic regions are particularly important loci because of their extreme positive correlation and pleiotropic effects, and could be considered for the identification of potential candidate genes corresponding to fruit bitterness, sweetness, firmness, and fruit size/weight of goji berries.

In this study, the whole genome resequencing (WGS) approach was employed on the two parents and interspecific 305 F₁ individuals and re-analyzed using in-house databases along with the L. barbarum reference genome. WGS of plant genomes offers profound insights into genome evolution and enables the discovery of a vast array of DNA markers, including SNPs, InDels, copy number variations, and presence/absence variations in crops. By improving the accuracy, output, speed, and affordability of genome-wide genotyping, WGS-based methods can transform QTL mapping (Huang et al., 2009; Xu and Bai, 2015; Pan et al., 2022);. Building an ultra-dense linkage map with precise genotyping using polymorphic markers such as SNPs has become feasible with the advancement of NGS technologies and the release of the goji berry reference genome. Nevertheless, there is a report of a goji berry high-density genetic map containing 15,240 SNP markers based on resequencing, and identified leaf- and fruit-related trait loci (Zhao et al., 2021). Similar studies utilizing a large number of SNP markers in different crops have reported the use of WGS to generate high-density genetic maps (Li et al., 2018; Guo et al., 2023). This will offer useful resources for marker-assisted selection, map-based gene cloning, and candidate-gene discovery. Our research presents an incredibly detailed genetic map of the goji berry, which was derived from whole-genome resequencing of both F₁ individuals and parental lines. The WGS method applied in this study for constructing genetic maps is notable for its integration of genotyping, SNP validation, and discovery (Xu and Bai, 2015). In our previous study, SLAF-seq was used to identify 3,495 SNP markers that were subsequently used to build a genetic map of goji without a reference genome (Rehman et al., 2020). To identify recombination breakpoints, we sequenced the parental lines at 49-fold and 43-fold sequencing depths, and the average sequencing depth of the F₁ individuals was 17-fold, yielding high-quality clean reads. Comparatively, the average depths in the female (NN) and male (YG) parents and F₁ individuals were 34-fold, 18-fold, and 2.58-fold, respectively (Zhao et al., 2021). Unmatched recombination intervals from previous genotypic data (Rehman et al., 2020) were filtered, and an ultra-dense genetic map was built using bins derived from reliable genotypic data. Following genotyping analysis, 12 linkage groups were assigned to 3,058 recombination bins, corresponding to 35,331 polymorphic markers. A previous study identified 8,507 SNP markers to create a consensus genetic map for goji berries (Zhao et al., 2021). Our findings demonstrate the effectiveness of the WGS strategy for both ultra-dense linkage map building and marker discovery, which accurately reflect the features of genomic and genetic variation in goji berries. Bin genetic mapping offers several benefits, including the ability to overcome the constraints imposed by conventional mapping techniques and capture the entire recombination landscape of a mapping population by selecting the most relevant markers from a large pool of marker datasets (Guan et al., 2022). During map construction, consistent markers paired as recombinant bins decreased marker placement faults and genetic distance estimation errors, while avoiding duplicate loci. This technique has proven effective in producing ultra-dense genetic linkage maps across a range of species (Li et al., 2011; Chen et al., 2014; Ma et al., 2016; Qin et al., 2022). Recently, several studies have constructed genetic maps of goji berries using NGS-based methods (Gong et al., 2019; Zhao et al., 2019; Rehman et al., 2020; Zhao et al., 2021; Yue et al., 2022). Gong et al. (2019) published a high-density genetic linkage map of goji berries using the double-digest restriction site-associated DNA (ddRAD) of genotyping by sequencing (GBS) method. The map contained 2,764 bins spanning a total length of 964.03 cM, with an average distance of 0.038 cM between adjacent bins. A high-density genetic map with a total genetic distance of 1702.45 cM was created using large-scale SLAF-seq for SNPs detection (Zhao et al., 2019). Using whole genome resequencing, Zhao et al. (2021) constructed an SNP-based consensus genetic linkage map with 12, linkage groups, 15,240 SNP markers, and a total distance spanning 3,058.19 cM with an average marker distance of 0.21. Earlier, a high-density genetic map of goji berry was developed using SLAF-seq, which covered a total distance of 1,649.03 cM and comprised of 3,495 SNP markers with an average interval of 0.47 cM between adjacent markers (Rehman et al., 2020). In the present study, a similar mapping population (Rehman et al., 2020) was employed, and a greater frequency of polymorphisms for map construction and notable variances in a number of features were caused by the large genetic diversity between the two parental lines. Using WGS, 3,058 bin markers including 35,331 polymorphic markers were identified and accommodated into 12 chromosomes The linkage map was 966.42 cm in total length, with an average distance of 0.03 cm between adjacent bins. Our genetic map was populated by a large number of additional SNP markers. In contrast, recent studies have shown that genetic maps of goji berries include fewer SNP markers (Gong et al., 2019; Zhao et al., 2019; Rehman et al., 2020; Zhao et al., 2021). Conversely, the resequencing-based linkage map revealed a smaller comprehensive map distance and a greater total number of SNPs, implying that WGS exhibits higher resolution than reduced-representation sequencing (Rehman et al., 2020). To date, an extremely dense and highly saturated linkage map of Lycium has been developed successfully. The collinearity of the genetic and physical maps was uniform for each linkage group, which is consistent with previous reports (Rehman et al., 2020; Zhao et al., 2021). Furthermore, visual assessment of the haplotype and heat maps of the genetic map, along with synteny analysis of the genetic map and reference genome, indicated that our map was produced appropriately and that the F₁ individuals were eligible for genetic study. This ultra-dense linkage map serves as a strong source for map-based cloning, marker-assisted selection breeding, and fine-mapping of potential genes for significant agronomic and economic traits.

In light of the significance of the important fruit quality and size/weight attributes of goji berries, phenotypic data of the studied traits were collected consecutively for three or four years. When evaluating data from different years corresponding to the eight traits, correlation analysis revealed a highly negative association between human-evaluated BT and ST. Conversely, a highly positive association was observed between ST and SSC, but no association was detected between ST22 and ST20/ST21, possibly due to the rejection of the null hypothesis for ST22 under normality tests ( Supplementary Table 3 ). In addition, fruit firmness was positively associated with fruit weight and length, whereas fruit size-related traits showed a highly positive correlation between different years, except for FD and FSI, with highly negative associations. The extremely negative correlation corresponding to fruit bitterness under different years with fruit sweetness, and SSC with FSI and FW, and highly positive correlation of fruit size-related traits might indicate a tight association among linked markers or even candidate genes due to pleiotropic effects. This information can help to determine candidate gene predictions. The outcomes of the correlation analysis between bitterness, sweetness, and SSC revealed a noticeable disparity between human taste evaluations and the readings obtained from digital instruments, such as refractometers. We posited that human sensory tests for fresh goji berries would exhibit lower precision and accuracy in determining bitterness and sweetness, potentially as a result of variations in taste buds and environmental factors, such as differences in planting areas. However, it is crucial to consider the sugar content determined by a refractometer when evaluating the taste of fresh goji berries. A similar study reported no significant association between color variations and visual judgment scores for bitterness and sweetness (Huang et al., 2021). In sensory evaluation, a negative correlation was found between the color score and bitterness, and a highly positive association was found between the color score and taste evaluation, indicating the influence of color differences on taste evaluation (Huang et al., 2021). The presence of flavonoid compounds with a dihydroquercetin structure is thought to contribute to the bitterness of goji fruit (Yuan et al., 2017). Moreover, the bitterness of wild peach fruits is associated with a higher concentration of phenolic compounds than that in landraces and improved cultivars (Predieri et al., 2005). The sugar-acid ratio and assessment of bitterness were linked to the titratable acid content in goji berry fruit (Huang et al., 2021). An intriguing taste pattern ranged from bitter to half-bitter, bitter after sweet, hemp bitter after sweet, sweet, slightly bitter after sweet, and whole bitter of freshly matured fruits was found in another study that assessed the germplasm of red goji berries (Yuan et al., 2022). The F₁ populations of Ningqi-7 and wild bitter goji berries displayed proportions of sweetness and bitterness, whereas the reverse cross produced only bitter F₁ individuals, indicating the heterozygosity of bitterness-related genes, which are controlled by dominant genes with additive genetic effects, whereas recessive genes may be responsible for controlling sweetness (Yuan et al., 2017). The distinct sweetness of L. barbarum goji fruit compared to L. chinense fruit is attributed to the presence of betaine, a naturally occurring amino acid. Research has demonstrated that higher concentrations of betaine are associated with increased sweetness (Lee et al., 2014).

In the present study, 99 genomic loci corresponding to eight physio-chemical and morphological traits were identified. The highest phenotypic variance explained by all detected QTLs varied from 1.15–16.95% with 2.52–29.25 LOD. The additive and dominance effects of each identified QTL were then quantified. The cumulative effect of an individual allele at a locus influences the trait value, which is known as the additive effect of a QTL. This is the result of the alleles assigned by each parent. The trait value is enhanced by alleles from the mother parent when there is a positive additive effect and enhanced by alleles from the father when there is a negative additive effect (Hu et al., 2021). In contrast, the dominance effect describes how one allele can supersede or mask the effects of another allele at the same locus; moreover, it also defines the interaction between alleles at a specific locus, which leads to deviation from an additive genetic model. Among breeding populations, the additive effect of a QTL is more significant as it directly contributes to the overall genetic variance than dominance effects; thus, it is considered a decisive factor in determining the genetic architecture of complex traits (Zhang et al., 2008). QTL mapping analysis found 17 total BT QTLs, two major loci with negative and positive additive effects and five overlapping genomic regions with the highest LOD scores of 29.25, and phenotypic variance of 16.95%. Similar studies have identified a bitterness QTL (qbt-c1–1) with a high LOD score in watermelons (Li et al., 2018) and 27 QTLs for fruit taste in peaches (Li et al., 2019). Two QTLs associated with bitterness (qBT.1, and qBT.2) were found in bottle gourds LG02 and LG09, respectively (Wu et al., 2019). Similar to earlier reports, goji fruit bitterness-associated loci were identified on chromosomes 1, 2, and 9 (qBT-1.1, qBT-2.1, qBT-9.1, and qBT-9.2) with a positive additive effect and the highest LOD score of 4.14 ( Supplementary Table 1 ). Twelve QTLs targeting ST, including a major-effect locus (qST-12.1) and five minor-effect QTLs, had the highest LOD of 14.12 ( Table 1 ; Supplementary Table 9 ). Similarly, 21 loci targeting SSC (°brix percentage), including qSSC-8.6 with the highest phenotypic variance (11.50%), and a positive additive effect were identified as the major locus ( Supplementary Table 9 ). Six other SSC loci overlapped with QTLs corresponding to FL and FF. An earlier study identified two fruit sweetness (SSC) QTLs (qFS3–1 and qFS5–2) with relative LODs of 5.48 and 9.2% PVE on LG03 and LG05, respectively (Rehman et al., 2020). However, qSSC-4.4 was identified (using the same 2019-year data) on chromosome 4, with a PVE of 4.14% and an LOD score of 2.86, exploiting an ultra-dense bin linkage map. Multiple SSC QTLs have been reported (Wu et al., 2021; Miller et al., 2022), including fruit flesh sugar content regulation QTLs (Qfru2–3, QBRX2–1, and QBrix6) on chromosomes 2 and 6 in watermelons (Sandlin et al., 2012; Ren et al., 2014). Consistent with previous findings, SSC QTLs (qSSC-2.3, qSSC-6.3, qSSC-6.4, and qSSC-6.5) were also detected on chromosomes 2, and 6 with a maximum LOD score of 27.24 ( Supplementary Table 9 ). Fresh fruit firmness is an important quality trait in goji berries that affects fruit texture, and is ultimately associated with shelf life and economic value. Fruit firmness is a quantitative trait regulated by multiple genes (Dai et al., 2022). Ten loci were associated with fruit firmness including a stable locus (qFF-7.8, qFF-7.13) ( Figure 5 ), and an overlapped loci with qFF-7.10 on chromosome 7, with LOD score of 3.01, and 5.90, respectively. Two fruit firmness QTLs (qff2.1, and qff5.1) were detected on chromosomes 2 and 5 in melons with high LOD scores of 3.8, 17.44, 28%, and 38% of variance, respectively (Dai et al., 2022; Chen et al., 2023). Similar investigations have found one locus (qff4.1) on chromosome 4 for sweet cherry fruit firmness (Cai et al., 2019) and blueberry fruit firmness QTLs on chromosome 8 (Cappai et al., 2020). One locus (qFF-5.2) corresponding to goji fruit firmness was identified on chromosome 5, with an LOD score of 3.02. Previously, a fruit firmness (FF) QTL (qFF10–1) was positioned on LG10 with a relative LOD of 4.81 and 8.1% phenotypic variance (Rehman et al., 2020). However, four FF loci (qFF-4.4, qFF-7.8, qFF-7.9, and qFF-12.5) were identified (using the same 2019-year data), and positioned on chromosomes 4, 7, and 12 with variance ranging from 2.04–7.08% and LOD scores of 2.53–7.53, respectively. It has been learned that the colocalization of QTLs governing various fruit size/weight-related traits provides an unusual opportunity to gain profound insights into the intricate mechanisms underlying fruit size traits. This colocalization results from the pleiotropic effects of particular genes, linkage disequilibrium, or convergence of several important genes within the same genomic area (Kato et al., 2000; Baytar et al., 2021). Moreover, because correlated traits are likely governed by comparable genetic causes, mapping them to similar genomic sites is expected to establish an interesting genetic architecture (Kato et al., 2000). In the present study, major-effect and colocalized QTLs targeting fruit size/weight accounted for a maximum LOD score of 11.54, and phenotypic variance of 12.23% was found to have positive and negative additive effects. Moreover, stable loci (qFW-6.6, qFW-6.9, and qFL-6.13) displaying negative additive effects were identified across the two environments with similar confidence intervals ( Table 1 ). Consistent with our findings, previous studies have reported major fruit weight QTLs (qAFW6.1, qAFW12.1) with positive and negative additive effects (Lin et al., 2010; Aruna et al., 2023). Nine fruit weight (FW) QTLs were positioned on LG10 with LOD ranging from 4.72–7.8 and 6.9–11.1% PVE in two individual years (2018–19) (Rehman et al., 2020). However, in the current study, 11 FW loci were identified based on the previous two years (2018–19) data and were positioned on different chromosomes with PVE ranged from 2.36–10.39% and LOD score ranging from 2.58–11.54 ( Supplementary Table 9 ). Similarly other goji studies have reported seven fruit weight loci in different linkage groups with a variance of less than 10%, including one stable FW loci on LG10 with a PVE of 59% (Zhao et al., 2019, 2021). Several loci corresponding to fruit length, along with two major- and co-localized loci (qFL-8.7, qFD-8.9; qFL-8.8, qFD-8.10) showing positive and negative additive effects with a PVE > 10%, were identified on different chromosomes ( Table 1 ; Supplementary Table 9 ). Previously, five fruit length (FL) QTLs were positioned on LG10 with LOD ranging from 6.08–8.1 and 8.8–11.5% variance, in two individual years (2018–19) (Rehman et al., 2020); nevertheless, 12 FL loci were identified and positioned on different chromosomes with phenotypic variance ranging from 1.15–11.16% and LOD score ranging from 2.53–22.70 ( Supplementary Table 9 ). A similar study on goji found six FL loci on LG10 and LG11 including two stable FL loci (qFL10 and qFL12), which explained up to 36% of the phenotypic variance (Zhao et al., 2019, 2021). Another investigation discovered two major fruit length QTLs on chromosomes 1 and 7 with a PVE of > 10% (Aruna et al., 2023). For fruit diameter Eight QTLs were detected, including a major and colocalized locus (qFD-8.10, qFL-8.8) with a phenotypic variance of 12.25%, revealing a positive additive effect ( Table 1 ; Supplementary Table 9 ). Earlier, QTL mapping results into five fruit diameter (FD) QTLs were positioned on LG10 with LOD ranging from 4.5–7.8 and 6.6–11.1% PVE in two individual years (2018–19) (Rehman et al., 2020); however, seven FD loci were identified and positioned on different chromosomes with variance ranging from 2.94–12.25% and LOD scores ranging from 2.59–17.08 ( Supplementary Table 9 ). Similar findings have identified four fruit width QTLs in different linkage groups, with PVE exceeding 10% for two QTLs (Zhao et al., 2019). QTLs targeting the fruit shape index were identified, with one major QTL (qFSI-5.8, qFSI-5.11) and two stable QTLs (qFSI-5.7, qFSI-5.10; qFSI-5.8, qFSI-5.11) harboring LOD values of 8.41, 18.99, and PVE of 6.07% and 14.02%, respectively. Previously, eight fruit shape index (FSI) QTLs were positioned on different linkage groups with LOD ranged from 4.55–7.08 and 6.6–10.2% PVE in two individual years (2018–19) (Rehman et al., 2020). Comparatively, four FSI loci were identified, which were positioned on chromosomes 2 and 5, with PVE ranging from 2.10–6.38% and LOD score ranged from 3.27–9.71 ( Supplementary Table 9 ). A similar study found 24 loci targeting the fruit index (FI) on different linkage groups with a PVE greater than 10%, and two stable QTLs (qFI7–1 and qFI7–2) corresponding to FI were detected in LG07, with a PVE of up to 19% (Zhao et al., 2019, 2021). The colocalization of fruit size QTLs, such as qFW-1.5, qFD-1.8, qFW-3.6, qFD-3.7; qFW-6.6, qFW-6.9, qFL-6.13; qFL-8.7, qFD-8.9; qFL-8.8, qFD-8.10; and qFL-6.14, qFD-6.15), may be associated with strong positive correlations among fruit weight, length, and diameter. The loci causing FW, FL, and FD were located at the same positions on chromosomes 1, 3, 6, and 8, supporting the high association between these three fruit size traits, suggesting that multiple genetic factors govern each of the three. The intriguing relationship between multiple loci that govern fruit bitterness, sweetness, firmness, SSC, and size/weight QTLs (qBT-8.2, qBT-8.3, qBT-10.1, qST-10.3, qBT-10.2, qFL-10.4, qBT-6.1, qFW-6.8; qFL-5.5, qBT-5.1; qSSC-4.3, qFL-4.5; qSSC-6.5, qFL-6.10; qSSC-12.3, qSSC-12.4; qSSC-7.5, 7.6, qFF-7.8, 7.10, 7.11; qFF-12.5, qFW-12.7; and qFL-2.3, qFSI-2.4) located on the same chromosomes might indicate a positive correlation, suggesting a functional interrelationship between these genetic regions ( Figure 5 ; Supplementary Table 9 ). This could provide a valuable resource for identifying candidate genes related to fruit quality and size in future studies. A similar phenomenon of colocalization and overlapping loci associated with different traits has been reported previously and is consistent with our findings (Xu et al., 2017; Aruna et al., 2023).