Identification of fruit size associated quantitative trait loci featuring SLAF based high-density linkage map of goji berry (Lycium spp.)

Background Goji (Lycium spp., 2n = 24) is a fruit bearing woody plant popular as a superfood for extensive medicinal and nutritional advantages. Fruit size associated attributes are important for evaluating small-fruited goji berry and plant architecture. The domestication traits are regulated quantitatively in crop plants but few studies have attempted on genomic regions corresponding to fruit traits. Results In this study, we established high-resolution map using specific locus amplified fragment (SLAF) sequencing for de novo SNPs detection based on 305 F1 individuals derived from L. chinense and L. barbarum and performed quantitative trait loci (QTL) analysis of fruit size related traits in goji berry. The genetic map contained 3495 SLAF markers on 12 LGs, spanning 1649.03 cM with 0.47 cM average interval. Female and male parents and F1 individuals` sequencing depth was 111.85-fold and 168.72-fold and 35.80-fold, respectively. The phenotype data were collected for 2 successive years (2018–2019); however, two-year mean data were combined in an extra year (1819). Total 117 QTLs were detected corresponding to multiple traits, of which 78 QTLs in 2 individual years and 36 QTLs in extra year. Six Promising QTLs (qFW10–6.1, qFL10–2.1, qLL10–2.1, qLD10–2.1, qLD12–4.1, qLA10–2.1) were discovered influencing fruit weight, fruit length and leaf related attributes covering an interval ranged from 27.32–71.59 cM on LG10 with peak LOD of 10.48 and 14.6% PVE. Three QTLs targeting fruit sweetness (qFS3–1, qFS5–2) and fruit firmness (qFF10–1) were also identified. Strikingly, various traits QTLs were overlapped on LG10, in particular, qFL10–2.1 was co-located with qLL10–2.1, qLD10–2.1 and qLA10–2.1 among stable QTLs, harbored tightly linked markers, while qLL10–1 was one major QTL with 14.21 highest LOD and 19.3% variance. As LG10 harbored important traits QTLs, we might speculate that it could be hotspot region regulating fruit size and plant architectures. Conclusions This report highlights the extremely saturated linkage map using SLAF-seq and novel loci contributing fruit size-related attributes in goji berry. Our results will shed light on domestication traits and further strengthen molecular and genetic underpinnings of goji berry; moreover, these findings would better facilitate to assemble the reference genome, determining potential candidate genes and marker-assisted breeding.


Background
Lycium barbarum L. and L. chinense Mill., two closely related Solanaceae species commonly known as goji, wolfberry or boxthorn. These are deciduous woody perennial shrubs, growing up to 4 m in height, widely distributed under warm and subtropical climate in East Asia, Southeast Asia and Europe [1][2][3]. It is used both in fresh and dried forms, while dry root bark and fruit also have functional importance in the medicinal industry particularly Traditional Chinese Medicine (TCM) [2]. Goji has been appreciated for its numerous health benefits as it increases the eye vision, inhibits the growth of cancer cells [4,5], anti-fatigue, anti-aging, enhance metabolism and ultimately improves the immune system [1,6]. Physiologically, the fruit of these plants contains many compounds of functional characteristics such as vitamin C, carotenoids, flavonoids and polysaccharides [3]. Ningxia as "Daodi" region of goji production, widely renowned in China. Somehow, due to recent market demand goji cultivation areas have stretched to new regions over different climatic zones covering from 82°E and 115°E to 30 o N and 45 o N [7]. Particularly, cultivation of L. barbarum L. is largely extended to semi-arid climate temperate continental; Ningxia, Gansu, Inner Mongolia, Plateau continental; Qinghai and Continental arid climate; Xinjiang [8], whereas L. chinense Mill. is notably confined to temperate monsoon climate like Hubei [7,8].
It is demonstrated that fruit weight controlling loci; fw2.2, fas and lc (WUSCHEL), contribute fruit size regulation by increasing carpel number, which influences yield acreage dramatically in tomato and other Solanaceae crops during domestication [9][10][11][12]. Also, fw2.2 ortholog in maize as cell number regulator gene (CNR) potentially increases heterosis vigor and yield [13]. It is witnessed that fruit size traits substantially contribute to yield components [10,13,14]. Previous genetic mapping and QTL studies in interspecific tomato population revealed fruit weight loci fw1.1, fw2.2, fw3.1, fw3.2 and fw11.3 each contributing more than 20% phenotypic variation. Some of these effective QTLs, fw2.2 [10], fw3.2 [15], and fw11.3 [16], responsible for fruit weight variations have already been cloned in individuals' studies. As shown, tomato fruit shape involved key QTLs, for example, Fasciated (fas) and lc (WUSCHEL) two important quantitative trait loci contributed synergistically to the locule number enhancement and collectively triggered fruit size in tomato [17,18]. Comparatively, L. barbarum fruit is highly valued due to its bright red color, sweet taste and high economic returns [7]. The fruit is 2 or 3 chambered, fleshy, juicy with few or up to 15 seeds per fruit. The fruit sizes differ significantly between L. barbarum and L. chinense; ZKLC1 fruit is approximately 2.5 cm in size than L. chinense of 1.2 cm. Current germplasm and varieties of the goji berry are diverse in their agronomic characteristics, especially oblong or narrow leaves and fruit size round to oblong [20]. We have learned from previous similar studies within the Solanaceae family that tomato [21], pepper [22], potato [23], and eggplant [12], have contributed extensively for conventional and molecular genetic research in terms of the development of molecular markers based high-density linkage maps, and mapping of qualitative and quantitative attributes [12]. Goji, a closely-related Lycium species in the Solanaceae, have been utilized for conventional genetic approaches, including hybrid development, studies on qualitative and quantitative traits [24]. Thus far, goji berry has been less explored in terms of fruit genetic aspects and breeding of quality parameters except for one recent report [25]. The accessibility of the high-density genetic map, established using genomic DNA, and single nucleotide polymorphism (SNP) markers, explore numerous opportunities to the QTL mapping, gene tagging, and isolation of candidate loci in goji berry as well as comparative mapping studies with other Solanaceae species [12,26]. Moreover, the presence of orthologous loci or genes have already been established in several reports [12,22,27]. Thus, we suggest that goji berry could explain the similar trend of genetic identities or represent a novel set of candidate genes triggering fruit size for future map-based cloning. However, larger fruit sized cultivars in goji berry industry is a major breeding goal and will be a breakthrough to further marketing expansion and production increase. Molecular breeding and markerassisted selection (MAS) approaches of goji berry have not yet been attempted; moreover, it is a perennial tree plant with long duration from germination to flowering.
Traditional breeding strategies limited due to intensive labor, cost and time inefficiency, which impeded cultivar improvement of goji berry [25].
Under rapid advancement of next-generation sequencing (NGS) technology, high-resolution genetic mapping and single nucleotide polymorphism (SNP) markers provide affordable tools for QTL mapping and MAS [28]. In the perennial fruit crops, the genetic map based on F 1 population provides a robust source to locate linkage between commercial traits and DNA markers [29]. Highly heterozygous mapping population can be successfully developed by interspecific crosses, and F 1 individuals determined under pseudo-testcross approach are widely applicable in perennial trees and forest herbs [30][31][32]. Recently, Gong et al. (2019) [33], constructed the first high-density genetic map of L. barbarum based on intraspecific F 1 population using ddRAD-seq, which contained 23,967 SNPs with the total genetic length is 964.03 cM and average interval is 0.040 cM. QTL analysis identified 8 significant loci targeting photosynthetic traits in 2 years dataset. Importantly, QTLs linked to net photosynthetic rate (P N ) and stomatal conductance proposed critical for the plant growth and development of goji berry [33]. SLAF-seq, a high throughput reduced representation sequencing strategy for large scale SNP discovery has been practiced due to numerous advantages [34,35], in different crop plants including pepper [36], maize [37], walnut [38], cauliflower [39], eggplant [40], and Salvia miltiorrhiza [41]. Moreover, SLAF-seq is widely applied to different woody perennial tree plants and efficiently generated genetic maps based on F 1 population [29,33,42]. Zhao et al. (2019) [25], reported wolfberry SNP based genetic map by SLAF-seq, which contained 6733 SNPs with the total length is 1702.45 cM and average inter-marker distance is 0.31 cM. Moreover, QTL analysis showed 55 QTLs, of which 18 QTLs detected for fruit index on LG11 and 2 significant QTLs for leaf index on two different linkage groups based on 2-and 3-years data, respectively [25].
Therefore, the current study executed SLAF-seq technology by establishing high-resolution linkage map utilizing 305 F 1 individuals of interspecific population derived from L. chinense Mill. cv. Daye and L. barbarum L. cv. ZKLC1. The map contained 3495 high-quality SLAF markers including 15,815 SNPs with a final distance of 1649.03 cM and 0.47 cM average marker interval. The average sequencing depth of the genetic map was 111.85-fold in the female parent, 168.72-fold in the male parent and 35.80-fold in F 1 individuals, accounting for greater depth than previously reported [25], which ensures the accuracy of markers along the genetic map. Notably, the QTL mapping analysis identified 42 significant genomic regions or potential molecular markers closely linked with corresponding fruit size-related traits.
The high-density linkage map and promising QTLs pinpointed in this study are prerequisite for marker-assisted breeding and provide possible underpinnings for isolation of potential candidate genes underlying fruit and leaves associated commercial traits in goji berry.

Morphological traits variability analysis
A total of 11 morphological traits, including fruit weight (FW), fruit length (FL), fruit diameter (FD), fruit shape index (FSI), number of fruits per end cluster (Nof/ec), number of fruits per node (Nof/n), number of seeds per fruit (Nos/f), 100 seed weight (100SW), leaf length (LL), leaf diameter (LD) and leaf area (LA), were evaluated based on 2 individual years (2018-2019) and extra year (1819). The analysis of variance (ANOVA) of morphological traits of 2 individual years showed significant differences (P < 0.05 or P < 0.01), which indicated variation among 305 F 1 progenies under different years (Table S1; Additional file 4). The descriptive statistical analysis showed coefficient of variation (CV%) ranged between 9 to 45 for FSI and Nof/ec in 2018, and 8 to 54 for FSI, FD and Nof/ec in 2019. However, 1819 data showed that CV% varied from 7 to 39 for FD and Nof/ec (Table S2; Additional file 5). The frequency distribution histogram and box chart among all morphological traits for individual years (2018-2019) and extra year (1819) showed normal distribution (Fig. S1a-z; Additional file 1). Normality test was performed using Kolmogorov-Smirnov (K-S) goodness of fit based on the absolute distance between cumulative distribution and values ranged from 0.02 to 0.11, suggested positive normal distribution of all traits evaluated. The correlation analysis (P < 0.05 or P < 0.01) revealed extreme significant positive association detected in a few comparisons among FL, FD, 100SW, Nos/f and FW; FD, FSI, 100SW, Nos/f, LL, LD, LA and FL; 100SW, Nos/f, LD, LA and FD; LL and FSI; LD, LA and 100SW; Nof/n, LL, LD, LA and Nos/f; Nof/ec, LL and Nof/n; LD, LA and LL; LA and LD or significant positive correlation between LD and FW; LL and FD; LD, LA and FSI; Nof/ec and Nos/f; LD and Nof/n. While highly significant negative correlation was detected between FSI and FD; Nof/n and 100SW; Nof/ec and 100SW ( Table 1). The extremely positive correlated traits might indicate tight association among linked markers or even candidate genes due to pleiotropic effects, such as FW and FL, FD, Nos/f, and FL and FD, FSI, LL, LD, LA. This information could help determine candidate gene prediction.
Analysis of SLAF-seq data generated from F 1 population The sequencing data were examined to ensure the validity of the SLAF library construction. According to the goji genome size and GC content information [43], the pepper genome was selected as the reference genome utilizing self-developed predicted software and established protocol [44]. For optimal enzyme digestion the following guidelines were practiced such as; 1) lowest presence of selected digested fragments in the repetitive sequences, 2) regular distribution of digested fragments in the reference genome, 3) length and specificity of the digested fragment and specific pilot experiment system must comply uniformity, 4) the final number of digested fragments or SLAF tags must meet the expected number of tags [44,45]. RsaI and HinCII were used for the SLAF library construction and 221,608 SLAF tags were obtained based on the pre-SLAF experiment.
The cross-check revealed normality of the SLAF library with paired-end mapping reads of 92.16% and enzymatic digestion efficiency of 93.11%. A total of 3021.32 M clean reads of data with a length of 200 bp were generated for goji berry, comprising 37 M, 30 M and 9 M total paired-end reads in the female, male parents and offspring, respectively. The source data quality was ensured by average Q30 ratio of 95.04 and 39.62% GC. Correspondingly, in the female and male parent total number of SLAFs developed were 155,002 and 217, 485 with an average depth of 72.43-fold and 60.43-fold of each SLAF marker, while for 305 F 1 individuals 199, 560 SLAFs were generated, with an average depth of 15.23-fold (Table 2). Of 494,472 high-quality SLAFs, 214,961 were polymorphic with 43.47% of polymorphism rate; furthermore, 279,356 non-polymorphic and 155 repetitive SLAFs were identified. After filtering out polymorphic SLAFs, a total of 40,616 SLAFs was classified into eight segregation patterns (Fig. 1). For the F 1 crosspollinated population, only 7 segregation patterns (ab×cd, ef × eg, hk × hk, lm × ll, nn × np, ab×cc, cc × ab) were retained for genetic mapping. We followed the genotyping quality criteria to filter out low quality and redundant markers using salient steps of filtering process (See "Methods" section). Collectively, we obtained 5669 total SLAFs out of which only 3495 high-quality SLAF markers were eventually utilized to construct an integrated (combined female and male parent map) high-density genetic map of goji berry with an average parental sequencing depth of more than 100-fold and offspring > 10-fold (Table 3). Among 3495 SLAF markers, 'SNP_only' were pre-dominant with 61.08% on the integrated map followed by 'InDel_only' (0.59%) and SNP & InDel (38.55%) (Fig. 2).

Basic characteristics of the genetic map
After performing fine filtering process out of 5669 SLAFs only 3495 high-quality SLAF markers that met the quality standard were mapped onto the integrated   [25], constructed the integrated genetic map based on intraspecific population with 23,967 SNP markers, spanned 964.03 cM of final genetic distance and 0.040 cM average marker interval, which reflected highly saturated map [33]. Similarly, another study reported wolfberry highdensity SNP based genetic map with a final distance of 1702.45 cM and 0.31 cM average interval [25]. In comparison with previous reports [25,33], the current study revealed interspecific genetic map of goji berry with final genetic coverage of 1649.03 cM and 0.47 cM mean marker interval. In particular, LG6 & LG10 was observed with least inter marker distances indicated maximum saturation and might be considered recombination hotspots in this population.

Assessment of high-density genetic map of goji berry
To evaluate the quality of the genetic map, we performed several mapping approaches; Firstly, the integrity of all mapped markers was carried out between 305 F 1  ). In the current study, the average percentage of the missing fragments were 0.0125%, represent the true quality of the map. It was verified that the interspecific F 1 population was highly purified and suitable for high-density genetic map construction. DNA markers located on each linkage group were dispensed regularly with an average interval of 0.47 cM between adjacent markers despite high-rate of recombination events in F 1 individuals. Thirdly, we displayed heatmaps to show the recombination frequencies among markers located on each linkage group to further evaluate the quality of the genetic map using pairwise recombination rate. The heatmaps depicted markers on the map were ordered accurately and pair-wise recombination rates were considerably low among adjacent markers and represented with yellow color diagonally. The purple color, particularly more visible in LG6 and LG10, reflected the highest recombination rates among adjacent markers ( Fig. S4a-l; Additional file 2). Finally, 18 markers with segregation distortion (P < 0.001) were integrated in the construction of the genetic map, accounting for 0.34% of all mapped markers.

QTL mapping analysis
By using the high-density genetic map of goji berry, we plotted a large number of QTLs covering 13 traits. QTL analysis was performed using the composite interval mapping model with MapQTL v. 6 Table 5). Other traits QTLs such as FSI, Nof/ec, Nof/n, 100SW, FS (SSC) & FF were located in different linkage groups ( Fig. S5, S6, S7). In addition, qFS-1, qFS-2 corresponding to fruit sweetness with maximum phenotypic variance of 9.2% and qFF-1 to fruit firmness with 8.1% of PVE were mapped on LG3, 5 and 10, respectively ( Fig. 4 and S6, Table 5). Collectively, 117 QTLs targeting 13 important traits were found (  (Table 5 and S3; Additional file 6). QTLs detected corresponding to each trait in 2 individual years were described as follows: Fruit weight (fw): Total 9 fruit weight QTLs were detected on LG10 and 6, of which major QTL was qFW10-1 with explained phenotypic variance of 11.1%. One stable QTL (qFW10-6.1) was located on LG10 with genetic distance from 65.89 to 71.86 cM and peak LOD value up to 6.01 along with corresponding PVE of 8.7% (Fig. 4). Moreover, 4 QTLs controlling FW was mapped on LG10 under extra year (1819) and qFW10-2 was confirmed as effective loci with maximum of 11.6% PVE (Table 5 and S3, additional file 6). Fruit length (fl): 5 QTLs responsible for fruit length were found on LG10. The maximum phenotypic variance of 11.5%, corresponding to qFL10-1 considered as significant QTL. qFL10-2.1 targeting fruit length was observed as stable QTL on LG10, spanned an interval ranged from 44.64-56.88 cM with corresponding highest LOD of 6.54 and 9.5% phenotypic variance (Fig. 4). Another 4 FL QTLs were detected on LG10 in 1819 with highest PVE of 12.4% corresponding to qFL10-1 (Table 5 and S3, additional file 6).
In our results, qLL10-1 was regarded as a major QTL with the largest explained 19.3% of variance and significant LOD of 14.21 (Table 5). Intriguingly, we encountered an interesting phenomenon of QTL co-localization among fruit and leaf related traits. Specifically, qFL10-2.1 was overlapped with tightly linked markers corresponding to LL, LD & LA among stable QTLs and represented a highly significant positive correlation (Tables 1 and 5). The proportion of variance (PVE%) between promising QTLs was > 10%, which indicated major contribution of stable QTLs controlling corresponding important traits. Moreover, we had noticed a number of co-located QTLs such as fruit weight QTLs  (1819), which also indicated highly positive significant correlation as well (Fig. 4, Tables 1 and 5).

SLAF based high-density mapping
A high-density and high-quality genetic map of goji berry assembled in the present study based on F 1 population derived from 2 different species of Lycium using SLAF-seq strategy, demonstrated to be a highly costeffective technique for the discovery of a large number of SNPs and wide-scale genotyping [34]. SLAF-seq has been utilized in many different crop species, even those without reference genome [33,42,46]. Gong et al. (2019) [33], reported very first goji berry high-density genetic map using an intra-specific F 1 population under ddRAD-seq. While the current perpetual report on genetic mapping using SLAF-seq will further enhance the technology and improve understanding of molecular genetics of goji berry. Though, SLAF-seq approach has rarely been utilized in goji berry to construct a genetic map for QTL analysis. Comparatively, previous GBS based approaches such as RAD and ddRAD-seq have limitation due to technical and procedural complexity, cost unfeasibility and lack of pre-designed scheme and deep sequencing for ensuring accuracy and efficiency [46]. Subsequently, SLAF sequencing takes advantage of high-throughput, cost-effective, deep sequencing technique based on double Barcoding system accompanying a larger population and pre-designed strategy for large scale SNPs, InDel detection and genotyping [34].
In the first place, we contrived a scheme along the foundation of well-assembled reference genome sequence. Due to unavailability of goji genome, we explored goji genetic and DNA content detail [43], and pepper was selected as a reference genome to carry out the sequencing and ensure the uniformity, density and efficiency of markers production. Under the pilot  [34], the recommended quality scores should not be less than Q30 and 6-fold minimal sequencing depth for each individual in SLAF-seq defined pipeline. Our dataset contained 3021.32 Mb paired-end reads and achieved 95.04% Q30 score. The sequencing depth for parents was 66.43-fold and 15.23-fold for individuals, which fulfilled the quantity and quality criteria of markers needed to conduct high-density genetic mapping. At last, we abided by strict principles of filtering to get high-quality SLAF markers (See "Methods" section). Two hundred fourteen thousand nine hundred sixty-one polymorphic SLAFs were obtained with a polymorphism rate of 43.47%, higher than previously identified [25,47,48]. In the current study, 8 segregation types were received with 24, 329 maximum number of markers accounting for (aa×bb) homozygous group, which was larger than previous reports [25,33]. A total of 5669 high-quality SLAF markers were obtained following filtration process. All the demonstrated outstanding features represented markers accuracy, high-throughput, efficiency at a feasible cost, support further use of the SLAF sequencing technique for goji berry.

High-density genetic map of goji berry
After building a successful SLAF library and the generation of high-quality markers.  [25]. Goji has gained popularity due to economic importance and its historical use as TCM [50]. Thus, molecular breeding and genetics aspects of goji berry still need to be unfolded. In the current study, 305 F 1 individuals were selected for high-resolution strategy, though the large number of individuals could bring higher recombination events [25,51], and accelerate the accuracy of fine QTL positioning [52].
The reliability of the genetic map was determined by haplotype maps, heatmaps, maximum gaps, mean gaps<=5, integrity of individuals on the map and segregation distortion (P < 0.001). The haplotype maps and heat maps generated for each linkage group revealed SLAFs distributed in a regular manner on all 12 LGs and correctly ordered. Moreover, pair-wise recombination rates were found lower significantly among adjacent markers except LG10 and 6. The maximum gap and average gap<=5 were larger than previous study [25]. LG1, 4 and 7 were found with varied marker density and maximum gaps greater than > 10 cM. Likely, low marker density areas depict strongly homozygous regions of the goji genome, a similar trend observed by Ollitrault et al. (2012) [53]. In contrast, high-density marker areas depict the centromeric location on the chromosomes contributing larger physical distance leading to lesser genetic distance. Relatively, high-density areas may interact with the genome at some position due to interspecific heterozygosity as described by Lind-  [53]. The integrity of the individuals on the genetic map was 99.03%, which reflected true high-quality and consistency of the map. Meanwhile, only 18 (0.34%) segregation distortion (P < 0.001) markers were observed in LG11 and 12, significantly lower than previously reported [25,42]. The segregation distortion among markers might be ascribable to the distant relationship of the parents and gametic or zygotic selection [55].
The current updated high-quality and high-density genetic map were found accurately designed and extremely saturated than earlier studies [25,41]. Therefore, we claim that the linkage map is precise and exhibit high-throughput, which would provide insights on the molecular biology of economically important medicinal tree, plant and subsequently strong grounds for QTL mapping and candidate gene identification particularly of fruit size-related traits in goji berry. Further, it will ensure to be an essential source for comparative genomic studies [56], and aid in assembling the reference genome of goji berry.

QTL mapping analysis for fruit size-related traits
QTL mapping using high-density genetic maps can provide a fruitful estimation in specific quantitative traits and data mining [36]. For conducting QTL analysis, different software is available limited to the crossing models such as IciMapping [57], MapQTL [58], and so on. The current study is based on F 1 interspecific population and cross-pollinated (CP) model, which is different than other population used for mappings such as RILs, back cross and double haploids [52]. Thus, MapQTL with composite Interval mapping (IM) is definitely a suitable software for conducting QTL analysis in CP model, likewise it has been successfully utilized in different species, pepper [52], cucumber [49], and walnut [38]. Zhang et al. (2010) [59], successfully plotted QTL analysis for 8 traits in sweet cherry. In a recent report by Gong et al. (2019) [33], a total of 32 QTLs were found corresponding to six photosynthesis-related traits using MapQTL in L. barbarum L. Similarly, in another study the interspecific high-density genetic map of goji berry has been explored for six agronomically important traits utilizing MapQTL [25].
We employed MapQTL v. 6.0 and characterized a large set of QTLs controlling fruit and leaf-related traits and fruit quality parameters as well. Particularly, 6 promising QTLs observed on LG10 and 12, targeting FW, FL, LL, LD and LA under 2 individual years consistently. For the extra year (1819), 36 QTLs were detected corresponding to FW, FL, FD, FSI, Nof/ec, Nof/n, Nos/f, 100SW, LL, LD and LA ( Fig. 4 and S6, Table 5). QTLs detected consistently in 2 years (2018-2019) dataset were regarded as stable QTL. Among several stable QTLs, fruit weight QTL (qFW10-6.1) was detected on LG10 with LOD up to 6.01 and explained phenotypic variance of 8.7%. In contrast, no fruit weight stable QTL was detected in previous study [25]. The fruit weight QTLs have been reported in different Solanaceae species, for example, in tomato 28 FW QTLs were identified. Specifically, fw1.1, fw2.1, fw2.2, fw3.1, fw3.2 and fw11.3 major fruit weight QTLs mapped in almost four interspecific individual studies and accounted for fruit weight variation during tomato evolution [21]. Similarly, 5 fruit weight QTLs were detected in pepper, and fw3.2 was found with significant impact and larger phenotypic variance (R 2 = 0.12-0.15) [22]. Pereira et al. (2018) [47], investigated QTLs regarding fruit morphological traits in Cucumis melo L. and FWQU5.1 was detected as a major QTL on LG5 controlling fruit weight with maximum PVE of 28.3%. Another study found, 3 QTLs targeting fruit weight in eggplant on different LGs with maximum phenotypic variance of 44% corresponding to qFW-9.1 [12]. For fruit length, we detected promising QTL (qFL10-2.1) on LG10 with peak LOD of 6.54 and 9.5% phenotypic variance. Zhao [29,40,60].
Among the leaf related traits, qLL10-2.1 stable QTL corresponding to leaf length was spanned on LG10, with peak LOD value and PVE up to 10.48 and 14.6%, respectively. Other stable QTLs (qLD10-2.1 and qLD12-4.1) controlling leaf diameter located on LG10 and 12 with the highest LOD of 6.64 and 5.91, and 9.5 and 8.5% phenotypic variance. The stable QTL responsible for leaf area was qLA10-2.1 spanned on LG10 with corresponding peak LOD of 7.74 and explained variance of 11%.
Nonetheless, previously only two stable QTLs were reported regarding leaf index (qLI10-2 and qLI11-2) on LG10 and 11, harbored 4 and 20 markers, respectively; while, no stable QTLs controlling leaf length and leaf width were reported [25]. Zhao et al. (2019) [25], detected total 55 QTLs targeting fruit weight, fruit length, fruit index, leaf width and leaf index except leaf length based on more than 2 years phenotypic data. Most of the QTLs were spanned on LG11 indicated cluster region for various traits. Moreover, 18 significant QTLs corresponding to fruit index were detected on LG11 covering an interval of 73.49-90.94 cM and one major QTL (qFL11-3) targeting fruit length was also observed with highest PVE of 30%. Among leaf related traits, 2 stable QTLs (qLI10-2, qLI11-2) targeting leaf index were detected on two different linkage groups for the 3-year dataset. Comparatively, in the present study, we detected 78 QTLs corresponding to all traits under investigation for 2 individual years of dataset along with 36 QTLs observed for an extra year (1819) spanned an interval ranged from 3.65-111.76 cM. In addition, 6 stable QTLs (qFW10-3.1, qFL10-2.1, qLL10-2.1, qLD10-2.1, qLD12-4.1, qLA10-2.1) were observed on LG10 and LG12 targeting fruit size-associated attributes (Fig. 4, Table 5).
LG10 was observed loaded with significant QTLs represented as a hot spot region for important fruit and leaf related traits. Interestingly, we found a couple of overlapping QTLs with those detected by Zhao et al. (2019) [25] such as qFSI1-1 and qFI1-2 [25], located both on LG1 and spanned an interval ranged from 123.73-124.76 cM & 114.81-123.77 cM, respectively. Another QTL set (qNof/ec5-3.1, qNof/n5-1.1 and qLW5) were overlapped on LG5 and covered genetic distance ranged from 68.04-76.61 cM [25] (Table S3, additional file 6). We also encountered colocalization phenomenon of qFL10-2.1 with qLL10-2.1, qLD10-2.1 and qLA10-2.1 along with tightly linked markers among stable QTLs (Fig. 4, S8a-o; Additional file 3). Additionally, it was observed that QTLs corresponding to FW, FL, FD, Nos/f, LL, LD and LA shared similar genetic position on LG10 throughout two individual years as well as in extra year. The co-localization of traits linked QTLs usually reflected positive significant correlation among different traits [22]. For example, fruit diameter and pericarp thickness analysis revealed 3 common positions on the genetic map of pepper (r = 0.86, 0.87), reflected dual genomic regions with common interaction for both traits [22]. So, we suggest that the current findings could reveal extra potential than previous similar studies and would provide a novel foundation to understand the genetics linked to fruit size-related attributes of goji berry.
It is one assumption that QTLs detected in the same genomic region targeting various traits may occur due to pleiotropy or linkage phenomenon [22]. Among promising QTLs, only 17 markers were found closely linked with FL, LL, LD and LA. It is also worth mentioning another set of 5 markers, which were located in stable QTLs related to FL, LL and LA on LG10. In the previous similar study, 3 markers were found significantly related to several traits in pepper, for example, 1 marker (LG8) influenced six traits, another marker (LG2) associated with four traits and third marker (LG3) related with five different fruit development traits [22]. Thus, we may hypothesize that those particular five genomic regions or loci may involve either a specific locus associated with pleiotropy, which might influence different aspects of goji fruit size and plant architecture or several linked loci indicates the specific influence on plant developmental responses. Further studies will require to target these promising QTLs for determining the likely candidate genes affecting different fruit development traits in goji berry. Furthermore, fruit sweetness (SSC) QTLs (qFS3-1, qFS5-2) were detected on LG3 and LG5 with relative LOD of 5.48 and 9.2% PVE and qFF10-1 responsible for fruit firmness were spanned on LG10 with peak LOD of 4.81 and PVE 8.1% (Fig. 4, Table 5). Ben Chaim et al. (2001) [22], reported similar events by identifying 3 QTLs in each year, but none as stable, corresponding to soluble solids concentration (SSC) in pepper with the largest phenotypic effect (R 2 = 0.12). Two QTLs (fi9.1, fi11.1) were detected controlling fruit firmness in pepper with 14% total variance [22], which indicated close agreement with our findings. Some other study conducted QTL analysis and found 2 QTLs (QBRX2-1, QBrix6) targeting flesh sweetness in watermelon [61].
However, transcriptome sequencing analysis is too indispensable to measure the expression pattern between different fruit and leaf developmental stages, which can exacerbate the mining process of candidate gene identification [25,49]. It has been brought out that promising orthologous QTLs or genes in similar crop species influence the same morphological or phenotypic impacts. Some physiologically associated trait comparison has been done practically and identified similar genetic positions [22,62]. Before the 90s, there was no QTL comparison has been made relating to fruit-producing crops until the identification of fruit associated traits in tomato [21,63]. Ben Chaim et al. (2001) [22], performed comparative QTL analysis between pepper and tomato and found quite interesting findings regarding fruit shape and fruit weight attributes. The orthologous QTLs in pepper and tomato targeting fruit weight on chromosomes 2 and 3 and one RFLP marker was found closely associated with fw2.1 and fw2.2 in pepper [22], and tomato [64], and another marker with fw3.2 and fw3.1 in pepper [22], and tomato [21], respectively. The above explanation could be helpful in the current study as Goji (Lycium spp.) being relative to the Solanaceae; thus, we might expect similar genetic identities or syntenic association, which existed in tomato and pepper [22,27], tomato and eggplant [12], and tomato and potato [23], which would be uncovered in later studies. Here, we have utilized common morphological traits related to the fruit and leaf development and these featured traits found to deliver a major contribution for domestication of fruit crops [12,56]. To further evaluate the potential of goji population and explore additional QTLs, we already have initiated developing recombinant inbred lines (RILs) by selecting full-sib mating with extremely distinguished offspring in 2019. RILs population have been successfully applied in many crop species such as pepper [52], peanut [65], and melon fruit [47]. Moreover, the enhanced QTL analysis could lead to exploring the genetic and molecular characterization of agronomically important fruit traits in goji berry.

Conclusions
Our findings of goji berry high-resolution linkage map using specific locus amplified fragment sequencing (SLAF-seq), may inculcate a powerful and more precise technique for de novo SNPs detection and high-density genetic map construction for species without reference genome sequence. A total of 3495 high quality SLAF markers was utilized for the construction of genetic map using 305 individuals from the F 1 interspecific population. The high-density genetic map provides 6 promising QTLs related to fruit weight, fruit length, leaf length, leaf diameter and leaf area. Remarkably, fruit weight, fruit length, fruit diameter, number of seeds per fruits and leaf related traits QTLs were shared similar location on the map, in particular, fruit length QTL (qFL10-2.1) were co-localized corresponding to leaf size-related stable QTLs with tightly linked markers on LG10. We might speculate that these promising QTLs have pleiotropic or linkage relationship controlling several traits. Our results will lay a solid ground for conducting positional cloning and marker-assisted breeding of goji berry. Moreover, our previous report results regarding the high-density genetic map of L. barbarum using intra-specific population and the current findings based on interspecific population may prove to be a better source for assembling the goji berry reference genome.

Plant materials
The planting material was comprised of Lycium chinense Mill. cv. "Daye" as the female parent and Lycium barbarum L. cv. "Zhongkelüchuan1 (ZKLC1)" as the male parent. Two cultivated varieties ZKLC1 and Daye were obtained from Goji germplasm reservoirs of Wuhan Botanical Garden, Chinese Academy of Sciences, Wuhan and South China Botanical Garden, Chinese Academy of Sciences, Guangzhou, China, respectively. The F 1 interspecific pseudo-testcross population was constructed by cross-fertilizing two parents in June 2016 and after sowing in the greenhouse during spring of 2017, the nursery was transplanted in the field at North-West China Bioagricultural Center, (38°28′05″ N, 106°16′23″ E) Yinchuan City, Ningxia Province, China. The female parent (Daye), also known as Chinese Goji, botanically has ovate to lanceolate leaves, shallow to dark purple flowers, fruit color range from dark red to orange, fruit shape oblong, fruit weight (0.38 g), fruit length (1.08 cm), fruit diameter (0.75 cm) (values are taken from collecting data), and seeds typically broader, rounded than L. barbarum with yellowish to yellow-brown color. The male parent (ZKLC1) has lanceolate to oblong leaves, bright to royal purple flower color, shiny red fruit color, round in shape, fruit weight (1.05 g), fruit length (1.79 cm), fruit diameter (1.28 cm), and seeds semispherical to flattish with light yellow color [66], (The parents Daye and ZKLC1 and progeny plant and fruit pictures can be seen in Fig. 5). In September 2018, 305 individuals randomly selected out of total 450 offspring were used for genetic linkage map construction and QTL mapping analysis, particularly of agronomic importance, fruit size and leaf related morphological traits such as fruit weight (FW), fruit length (FL), fruit diameter (FD), number of fruits per end cluster (Nof/ec), number of fruits per node (Nof/n), number of seeds per fruit (Nos/f), 100 seed weight (100SW), and leaf length (LL), leaf diameter (LD), leaf area (LA). Despite morphological traits, fruit quality traits such as fruit sweetness-soluble solid contents o brix FS (SSC) and fruit firmness (FF) were also included based on one-year dataset (2019).

Phenotypic traits estimation
The morphological trait data were gathered during two crop seasons, the summer (July) of 2018 and 2019. For fruit and leaf related traits: FW, FL and FD, LL, LD and LA, 20 fully ripened fruits and 10 fully grown first true leaves were picked randomly from progeny and parents' plant from the field, then photographed on scale paper (35 × 25 cm) by using Nikon D750 DSLR (Nikon Corporation, Japan) in the year 2018. The digital images were processed to measure the above-mentioned traits by ImageJ [48], software v. 1.52a [15,16]. Conversely, for 2019 summer season, a quicker yet precise method was adopted to assess the fruits and leaves related traits. Firstly, a digital scanner (MICROTEK ScanMaker i600) was employed to photograph the fruits and leaves, then built-in Chinese software Ban Ben, WSEEN (v 2.1.2.4) used to measure FL, FD, LL, LD and LA according to the developer's instructions. FSI was measured by dividing FL with FD. The other trait data, such as Nof/ec and Nof/n were collected visually in the field based on a maximum of 10 repeats per individual plant. Nos/f and 100SW were measured in the laboratory by taking 5 repeats and after overnight complete drying

DNA extraction and SLAF library construction and sequencing
Prior to DNA extraction, the young leaves from parents and 305 progenies were collected and dried with desiccant (silica gel). All the leaves DNA were extracted with Plant Genomic DNA Kit (TIANGEN BIOTECH (BEIJING), CHINA, CO., LTD.) according to the manufacturer's instructions with slight modifications. The quantity and quality of the extracted DNA samples were evaluated with a NanoDrop ND-2000 spectrophotometer (Thermo Fisher Scientific) and electrophoresis with 1% agarose gel, respectively. Specific length amplified fragment sequencing (SLAF-seq) was performed to genotype the 305 individuals from F 1 population and two parents according to Sun et al. (2013) [34], with minor modifications. The SLAF-seq strategy framework was adopted as the genomic DNA of each sample was subjected to restriction enzyme digestion using RsaI and HinCII (New England Biolabs, NEB, USA). The excised fragments were integrated with a single nucleotide A using Klenow Fragment (3′ → 5′ exo-, NEB), ligation of dual taglabelled sequencing adapters to the A-tailed fragment was performed with T4 DNA ligase (PAGE purified, Life Technologies, USA) [52]. The PCR amplification was completed by using diluted DNA samples after ligation, Q5®, HF-DNA polymerase (NEB), dNTPs, and primers set (forward primer: 5′-AATGATACGGCGACCA CCGA-3 and reverse primer: 5′-CAAGCAGAAGACGG CATACG-3). The target fragments (364-414 bp) were purified using a QIAquick gel extraction kit (Qiagen, Hilden, Germany). The purification of the PCR products was carried out using Agencourt AMPure XP beads (Beckman Coulter, High Wycombe, UK) and then pooled. The 2% agarose gel electrophoresis was used for pooled sample separation. After diluting gel-purified products, pair-end sequencing with a standard of 125 bp was performed using an Illumina Hi-Seq 2500 system (Illumina, Inc., San Diego, CA, USA) to sequence the SLAF in the quality-tested library at Biomarker Technologies Corporation, Beijing, China [44]. We used Japanese rice, Oryza sativa L. japonica genome as control with a genome size of 382 Mb (http://rapdb.dna.affrc.go. jp.), which claimed to verify the reliability and validity of the testing process and followed the same treatments in accordance with the goji mapping population.

SLAF sequencing data analysis and segregation pattern grouping
The grouping and genotyping analysis of SLAF sequencing data was performed in detail as described [34]. The raw reads were identified by sequencing to locate reads of the samples through dual index barcode sequences and low-quality reads of quality score < 20e were filtered out. The clean reads were clustered together through BLAT software [67], based on the similarity of more than 90%. After trimming of barcodes and the terminal 5-bp positions from each high-quality read, which later defined as one SLAF locus [66]. From each SLAF locus SNP loci were detected among parents, and SLAFs greater than 3 SNPs filtered out. According to the parents and offspring, reads with sequencing depth > 140.28-fold and > 35.80-fold were used to define alleles corresponding to each SLAF locus, respectively. One SLAF locus can hold mostly 4 genotypes for diploid species, hence more than 4 alleles in SLAF loci considered as repetitive SLAFs and discarded later. SLAFs with 2 to 4 alleles were recognized as polymorphic and known as potential markers. All SLAFs polymorphism loci were clustered into 8 segregation patterns as follows: (ab×cd, ef × eg, hk × hk, lm × ll, nn × np, aa×bb, ab×cc, and cc × ab). As the mapping population was F 1 interspecific from two highly heterozygous parents, thus all patterns were required to construct a high-density genetic map except (aa×bb).
A Bayesian approach was utilized for genotype scoring and quality as demonstrated [34]. To identify highquality SLAF markers 3 steps were taken into account. Firstly, the markers with average sequencing depths in each offspring found > 35.80-fold and parents > 140.28fold. Secondly, markers exceeding 10% missing data were discarded. Thirdly, markers with significant segregation distortion based on chi-square test (P < 0.001) were disqualified, but added later as additional markers during high-density genetic map construction.

Construction of high-density genetic linkage map and evaluation
Primarily, the marker loci were segregated into 12 linkage groups (LGs) under modified logarithm of odds (MLOD) scores > 5. To construct high-quality and highdensity map reliably, an advance high map strategy was required to assemble SLAF markers in a particular sequence and exact genotyping errors within linkage groups [19]. The maximum likelihood method by Van Ooijen [68], was employed to construct a high-density genetic map of goji berry and SMOOTH algorithm was utilized for genotyping error correction [69]. K-nearest neighbor algorithm was used to ascribe missing genotypes [70], and Kosambi mapping description was applied to estimate the genetic map distances between adjacent markers in centimorgan (cM) [71]. The genetic map was evaluated further following integrity of genetic maps, haplotype map, heat map and segregation distortion analysis (P < 0.05).

QTL mapping analysis
The QTL mapping analysis was executed for a mean of 11 fruits and leaf size-related traits and 2 fruit qualityrelated attributes and analyzed by MapQTL v. 6.0 [58], which employed inclusive composite interval mapping (ICIM) to detect QTL loci among 12 linkage groups for an integrated map of goji berry. A 1000 permutation test (PT) was performed to determine the logarithm of odds (LOD) for assessing statistical significance of each QTL [72]. For reliability of determining QTLs with major and minor effect, different LOD scores were assumed. At first, a LOD threshold corresponding to 0.99 confidence was considered, and if there was no mapping interval, the threshold value corresponding to 0.95 confidence was used; if there was still no positioning interval, then the LOD threshold value of 0.90 confidence was considered. Lastly, if there was still no mapping interval detected, then the permutation test result was not carried out, and the threshold value was manually lowered to 3.0, 2.5, and 2.0. The phenotype variance percentage as elucidated by its corresponding QTL (phenotypic explained variance %) was estimated within the segregating population on the basis of population variance using MapQTL. The QTL naming criteria were followed as described by McCouch et al. (1997) [73]. The QTLs positioning and genetic map were assembled by BioMercator v. 4.2.1 [74]., and MapChart v. 2.2 [75].

Statistical analysis
The general descriptive statistical analysis, including frequency distribution was carried out on Origin Pro v. 2016 (https://www.originlab.com/). The Analysis of variance for 2 individual years and Pearson correlation were accomplished on SPSS v. 25 software package (SPSS Inc., Chicago, IL, United States).
Additional file 1: Figure S1 Figure S2 The marker integrities of each individual in mapping population. *The x-axis indicates all 305 individuals along with the specific code name, while the y-axis shows markers integrity.
Additional file 2: Figure S3(a-l). Haplotype mapping of 305 F 1 individuals based on 12 linkage groups of integrated maps. *In each map, horizontal line stands a marker and column shows a chromosome in a sample. The first column of the map shows a paternal chromosome and the second column as maternal chromosome, and individuals are separated by blank columns. Green color indicates first allele from the parent, blue as the second allele from the parent, white color shows not judged event, and grey indicates a missing event. The position where color changes in the same column display reorganization events. Figure  S4(a-l). Heatmaps of 305 F 1 individuals based on 12 linkage groups of integrated maps. *Each row and column represent marker arranged in the orderliness of the map. Every small square show recombination rate between two markers. The color change trend from yellow to red to purple display reorganization rate from small to large. Yellow color indicates closer marker recombination rate, whereas closer the color from yellow to purple farther becomes the recombination rate.
Additional file 3: Figure S8  Additional file 7: Figure S5. QTL mapping on integrated linkage map in 2018. *The different color pattern represents identified QTLs for agronomic traits of goji berry such as FW_fruit weight, FL_fruit length, FSI_fruit shape index, Nof/ec_ number of fruits per end cluster, Nof/n_ number of fruits per node, Nos/f_ number of seeds per fruit, 100SW_100 seed weight, LL_leaf length, LD_leaf diameter, LA_leaf area.
Additional file 8: Figure S6. QTL mapping on integrated linkage map in 2019. *The different color pattern represents identified QTLs for agronomic traits of goji berry such as FW_fruit weight, FL_fruit length, FSI_fruit shape index, Nof/ec_ number of fruits per end cluster, Nof/n_ number of fruits per node, Nos/f_ number of seeds per fruit, 100SW_100 seed weight, LL_leaf length, LD_leaf diameter, LA_leaf area, FS_ fruit sweetness ( o brix) and FF_ fruit firmness Additional file 9: Figure S7. QTL mapping on integrated linkage map in 1819. *The different color pattern represents identified QTLs for agronomic traits of goji berry such as FW_fruit weight, FL_fruit length, FSI_fruit shape index, Nof/ec_ number of fruits per end cluster, Nof/n_ number of fruits per node, Nos/f_ number of seeds per fruit, 100SW_100 seed weight, LL_leaf length, LD_leaf diameter, LA_leaf area.