Dissecting the Genotypic and Environmental Factors Underpinning the Quantitative Trait Variation in a Set of Wild Tomato ( Solanum habrochaites LA1777) Introgression Lines

: Exotic libraries have proven a powerful tool for the exploitation of wild relatives and quantitative trait loci (QTLs) detection in crop species. In early 2000, an introgression line (IL) population of the wild tomato Solanum habrochaites (SH) acc. LA1777 was developed and made publicly available. Despite the potentiality of the donor parent, so far, these lines have been poorly explored for their agronomic performance and for the identiﬁcation of genomic regions underlying the variation of quantitative traits (QTLs). Here, we report the evaluation of 19 morpho-agronomic and chemical traits on a set of 39 ILs grown in three consecutive ﬁeld seasons with the aim to: (a) Determine the overall phenotypic performances of the studied collection, (b) estimate the inﬂuence of the genotype (G) and the year of cultivation (Y) and their interaction on the traits analyzed, (c) investigate the plasticity of the traits, and (d) identify whole-genome QTLs in the wild SH background. The ILs showed lower productivity compared to the control genotype, while no major effects were found for the morphological fruit-related traits. Instead, a general increase in the soluble solids content was observed. The combined analysis of G × Y highlighted a major effect of the genotype on trait variation, although yield-related traits were more inﬂuenced by environmental factors. In total, 75 associations for 17 traits were detected. Major QTLs increasing soluble solids, pericarp thickness, and trichome density were respectively found on chromosomes 1, 5, and 11 with a percentage variation (PV) of 24.01%, 32.49%, and 200%. Furthermore, different QTLs increasing the color intensity and fruit shape were detected. These results suggest that SH could be a potential source of favorable alleles for qualitative traits despite its inferior phenotype compared to the cultivated parent. The evaluated set of SH LA1777 ILs is a potential for novel allele discovery in wild tomatoes and for breeding purposes towards the exploitation of the available introgressions and for the pyramiding of traits.


Introduction
Tomatoes (Solanum lycopersicum, Solanaceae; 2n = 2x = 24) are considered one of the world's most important vegetable crops widely spread with a production estimation of 180 million tons and particularly appreciated for their nutritional properties [1]. As a model in plant research, the tomato has been pioneering for exploiting exotic germplasms toward the dissection of the genetic factors underpinning the variations of complex traits. Since the early nineties, the interspecific hybridization of S. lycopersicum with wild relatives has led to the development of various recombinant mapping populations able to enhance the identification of loci involved in agronomic and quality traits, as well as in biotic and abiotic stresses resistance [2,3]. Among the different types of experimental populations, introgression lines (ILs) were the most successful for allele transfer in precision breeding, providing, furthermore, an excellent material for quantitative trait loci (QTLs) mapping. ILs harbor a single marker-defined chromosome segment from the donor wild parent in the cultivated background [4]. The accuracy for mapping polygenes relies on the direct association of the phenotypic variation to the introgressed segment, minimizing the epistatic effects from other wild genome regions. Ideally, an IL population is supposed to cover the whole wild genome and is a permanent resource, being developed after subsequent backcross and self-fertilization cycles. This allows tests across multiple locations and seasons, providing a powerful tool to estimate the genotype (G) and the environmental (E) effects.
In fact, dissection of the G and E factors and their interaction is essential for determining those traits highly inherited in the offspring and for the estimation of phenotypic plasticity as the capacity of a genotype to express different phenotypes according to different environmental conditions [5]. Therefore, ILs play a key role for successful genetic improvement programs, facilitating, furthermore, the pyramid breeding of favorable genes [6]. In this context, the genetic variability enclosed in wild relatives can be exploited to address the future challenges of agriculture, such as the growth of the global population, the adaptation to climate changes, and the demand for quality and sustainability for food supply [7]. In recent years, ILs and near-isogenic lines (NILs) sets were developed in tomatoes through hybridization with different wild species, including S. pennellii [8], S. habrochaites [9][10][11][12], S. lycopersicoides [13,14], S. pimpinellifolium [15], and, more recently, S. sitiens [16]. Solanum habrochaites S. Knapp and D.M. Spooner (previously, Lycopersicum hirsutum Dunal) is a potential source of beneficial alleles in tomato breeding, including quality traits [17], chilling tolerance [18,19], and disease resistances [3,20]. There are over 90 accessions distributed in the native origin area for the species [21]; among these, acc. LA1777 is a green-fruited accession widely studied for resistance to the main tomato and potato pathogens, including begomoviruses and fungi [22][23][24]. It is also reported as resistant to insect pests, thanks to the presence of trichomes and of secondary metabolites produced by their glandules [25][26][27]. Finally, this wild species is a valuable source of beneficial alleles involved in fruit ripening [28] and an increase of soluble solids [10,29,30]. Monforte and Tanksley [9] developed a set of 99 near-isogenic and backcross recombinant inbred lines from a cross between a processing tomato cultivar E6203 and S. habrochaites LA1777. The population provides coverage of more than 85% of the wild genome based on genotyping with~93 RFLP (Restriction Fragment Length Polymorphism) markers. Despite that the population is publicly available, the enclosed potentialities for mapping plant and fruit traits have not yet been broadly investigated. Indeed, the identification of QTLs and genes affecting yield-related components have been, so far, mainly focused on the S. pennellii ILs, which have been extensively exploited for the identification of over 3000 QTLs [31,32] and for map-based cloning of the first two QTLs in tomatoes: fw2.2 for fruit weight [33] and brix9-2-5 for the soluble solids content [34]. Attempts toward the assessment of the agronomic performance of a subset of LA1777 ILs in hybrid combinations have been reported [35], although no QTL analysis has been performed. Instead, specific ILs were used for the fine mapping of QTLs for agronomic traits on the distal parts of chromosomes 1 and 4 [10,11]. Therefore, a lack of information at the whole-genome level still occurs for the LA1777 introgression library.
In the present study, we report the evaluation of 39 ILs for nineteen agronomic, morphological, and quality-related traits over three consecutive field seasons with the following objectives: (i) Determine the phenotypic performances of the selected set, (ii) estimate the genotypic and environmental effects and the plasticity of the traits, and (iii) identify wild introgressions responsible for the quantitative variations of the traits studied. The set of ILs was chosen based on the minimal number of introgressions able to cover, as much as possible, the whole genome, in agreement with previously reported data [9,36].

Plant Material
Plant material consisted of 39 S. habrochaites (acc. LA1777) (hereafter, SH) introgression lines and by the recurrent parent S. lycopersicum cv E6203 (hereafter, SL). The population was developed through an advanced backcross scheme derived from the cross between SL and SH. From the generation of BC 2 S 3 families and further marker-assisted selection, lines containing a single or few segments of the wild relative in the cultivated background were obtained [9].
The set of ILs was chosen in order to cover as much of the wild genome based on the published information [9]. Seeds were obtained by the C. M. Rick Tomato Genetics Resource Center (TGRC, University of California, Davis, CA, USA) [36] (Table 1).

Field Trials and Phenotypic Characterization
Plants were grown in independent trials across three consecutive seasons at the experimental farm of the Research Centre for Vegetable Crops located at Battipaglia (Salerno, Italy). Geographical and pedological characteristics of the site and climatic conditions of each season are reported in Table 2. Field tests were conducted using a completely randomized experimental design with 10 plants from each IL-SH and 40 plants from the control SL. Each year, the sowing was done in March, and the seedlings were transplanted in April in single rows, adopting distances of 100 cm between the rows and 30 cm along the rows. For plant fertilization, 200 kg ha −1 of N, 200 kg ha −1 of P 2 O 5 , and 120 kg ha −1 of K 2 O were applied according to what was established for the same site [37]. Irrigation scheduling was based on crop evapotranspiration. A total of 270, 260, and 290 mm of water were applied in 2014, 2015, and 2016, respectively, by drip irrigation, restoring 40% of the total available water depleted. The effective control of diseases and harmful insects was pursued, with the aim of having healthy plants until the end of the cycle. Other cultivation techniques included stakes as support and galvanized wires [38].
Agronomical and morphological traits were assessed when 95% of the plants of the control genotype (SL) had ripened fruits. Yield-related traits were performed on each plant and included: total yield (TY) (in kg), assessing the total weight of fruits (green and red), red yield (RY) (in kg) as the weight of whole, ripened fruits at commercial status, and green yield (GY) (in kg), including all unripe fruits at harvest time (in kg). Uniformity of fruit maturity (UM) was then measured as follows: The weight of the vegetative part (PW) (in kg) was obtained weighing the epigeal portion of the plants after harvesting all fruits; biomass (BM) (in kg) was measured as: The fruit traits assessment included the average weight (FW) (in gr) obtained by dividing the total yield by the number of fruits harvested, the length (FL) and the width (FD) (in cm) measured on 10 fruits/plant by using a ruler, and the fruit shape (FS) index as follows: where values < 1 indicate fruits more flattened, and values > 1 indicate fruits more elongated. External (EC) and internal (IC) colors were measured by scale (1 = yellow fruits and 5 dark red fruits), as well as pericarp thickness (PER) and puffiness (PUF) (1 = thin pericarp, low puffiness and 3 = thick pericarp, high puffiness). Chemical traits were measured on a bulk of 10 well-ripened fruits/plants. Fruits were washed and dried and then sliced and homogenized in a Waring blender (2-L capacity; Model HGB140, Waring Commercial, Stamford, CT, USA) for 1 min. Soluble solids content (SSC) was measured using a digital refractometer (Refracto 30PX, Mettler-Toledo, Novate Milanese, Italy), and the results were expressed as • Brix per 100 g of fresh weight (fw). The pH was determined using a pH-Matic 23 titroprocessor equipped with a pH electrode with a temperature sensor (model 5011T) (Crison Instruments, Barcelona, Spain).
The estimation of the yield of soluble solids produced per plant (BY) (in gr) was estimated as follows: Trichome density (TRIC) for each plant was estimated at the base of the stem and on the pedicel of flowers using a visual scale (1 = no or low density of trichomes and 3 = high density of trichomes).

Data Analyses
All phenotypic traits were subjected to an analysis of variance (ANOVA) test. Trait means for the control and the set of ILs were compared by using Tukey's HSD (honest significant difference) test. Results with p less than 0.05 were considered statistically significant. Broad-sense heritability (H 2 ) was estimated as: where "VG" is the genotypic component of variance, and "VP" is the total phenotypic variance. The coefficient of variation (CV) in percentage was calculated as: in which "SD" is the standard deviation, and "m" is the mean for a trait. A two-factorial linear model was used to determine the main effects of the genotype (G), season of cultivation (Y), and their interactions (G × Y) for each trait. Mean square values (MS) were used to estimate the magnitude of the observed effect, while the total sum of squares in percentage (TSS%) was calculated by dividing the TSS of the variable by the total TSS.
The coefficient of relative phenotypic plasticity (CRP) for each trait was calculated according to the formula of Dingemanse et al. [39] as follows: where "Vi" is the variance of the individual (i), and "VPp" is the overall phenotypic variance of the population. The CRP index close to 0 indicates no plasticity (high stability); on the contrary, higher values indicate an increase of plasticity (decrease of stability). Correlations across the genotypes for phenotypic traits were calculated using the Pearson's test at p < 0.01 after Bonferroni's correction for multiple comparisons. The correlogram was constructed and visualized using the Corrplot package implemented in R version 3.0.2 (R Development Core Team). Principal component analysis (PCA) was carried out to determine what are most effective descriptors in discriminating among accessions using the computer package XLSTAT 2012.1 and visualizing the similarities among accessions. Experimental data were statistically elaborated using R. The relationships between pairs of phenotypic and molecular data matrices were computed by the Mantel test using Pearson's r-value [40].
Chromosomal regions responsible for the quantitative variation for traits measured were detected comparing the mean phenotypic values of the tested ILs in each experiment to the recurrent parent SL using Dunnett's test with Type I error α ≤ 0.05 [41]. When an IL showed significantly different effects compared to the control in one or more seasons, putative quantitative trait loci (QTLs) were assumed in the chromosomal region covered by the IL. If a QTL was detected in multiple ILs, the likely position was inferred on overlapping chromosome segments. For each IL, the percentage variation (PV) compared to the control, and underlying a QTL region was calculated as: where, for each trait, "m" is the mean calculated across all replications.

Phenotypic Variation and Trait Performances of SL and SH ILs
In Table 3, the ANOVA results and the parameters mean, minimum, maximum, coefficient of variation (CV), and broad sense of heritability (H 2 ) specify the performance of the recurrent parent SL "E6203" and the 39 SH ILs set for each trait. Considering the average of three years of cultivation, no significant differences between the sets of ILs and SL were found for GY, UM, PW, IC, EC, and TRIC, while PUF and SSC, as well as FS, differed at p < 0.05 and p < 0.01, respectively. In average, the wild SH alleles conferred less productivity when introgressed in the SL background, while no major effect was noticed for the morphological fruit-related traits. Instead, a general increase in the soluble solids content was detected in the tested ILs, although no effect was found for the production of the soluble solid. Across all years, the CV% was lower than 36% for fruit traits in both SH ILs and the recurrent parent. A higher variation was shown by yield-related traits, reaching peaks for GY. Overall, the lowest coefficient of variation was found for FS, highlighting very low or absent morphological segregation within the ILs. The broad-sense heritability ranged from 0.84 (GY and CI) to 0.96 (FS). The fruits and chemical traits showed, on average, a greater heritability ( Table 3).
The mean, ranges, CV, and results of the post hoc Tukey's considering each independent season in both SL and the set of SH ILs are reported in Table S1. Highly significant differences (p < 0.001) between the ILs and the recurrent parent were found for 17 out of the 19 traits analyzed. Only the fruit shape and trichomes differed at p < 0.01 and p < 0.05, respectively (Table S1).
The average phenotypic values and the percentage variation compared to the control revealed a large variability of the studied ILs set for the considered traits (Tables S2 and S3).
In all trials, the introgression lines were less productive than the control. Overall, a decrease of the yield was observed in both E6203 and the ILs from the first to the third season. The same trend was observed for plant weight and biomass, with the lowest values exhibited during Y3. On the contrary, the weight of the fruits was higher during the second year of cultivation with respect to the other seasons, showing greater values in the control genotype across all trials (Table S2). The soluble solids content was found to be higher than the control in almost all lines during the three-year evaluation period ( Figure 1). Only for Y2, variable levels were observed. Overall, the soluble solids values reached contents up to +55% with respect to E6203.

Genotypic and Environmental Components Underlying Trait Variation
The results of the combined analysis of variance for the traits evaluated across three years of cultivation are given in Table 4. A high level of significance (p < 0.01) was found for the genotypic component (G), while the variation due to the cultivation season (Y) was lower for trichomes and fruit characteristics showing p < 0.05 for FS and CE, and there was no significance for the CI. As for the G × Y interaction, highly significant differences (p < 0.01) were found in all traits, except for FW and pH, which showed a p threshold < 0.05 and no significance, respectively.
For most traits, the main source of variability was due to G, which accounted, on average, for 61.90% of the total variation, expressed by TSS%, ranging from 12.35% for GY to 90.08% for TRIC. The traits under a strong genotypic effect were FL, CE, and CI, exhibiting TSS values above 80%. On average, G accounted for 43.01%, 75.28%, and 60.94% for yield-related traits, fruit characteristics, and chemical features, respectively. The partitioning of the TSS% in the other components affecting the variation indicated that Y and G × Y accounted for, on average, 16.47% and 21.77%, respectively, resulting in, generally, less influence of Y and G × Y than G on the analyzed parameters.
Among the assessed traits, GY and RY were the traits most affected by seasonal fluctuations, showing a TSS of 65.21% and 37.72%, respectively, while for all fruit traits, except

Genotypic and Environmental Components Underlying Trait Variation
The results of the combined analysis of variance for the traits evaluated across three years of cultivation are given in Table 4. A high level of significance (p < 0.01) was found for the genotypic component (G), while the variation due to the cultivation season (Y) was lower for trichomes and fruit characteristics showing p < 0.05 for FS and CE, and there was no significance for the CI. As for the G × Y interaction, highly significant differences (p < 0.01) were found in all traits, except for FW and pH, which showed a p threshold < 0.05 and no significance, respectively.
For most traits, the main source of variability was due to G, which accounted, on average, for 61.90% of the total variation, expressed by TSS%, ranging from 12.35% for GY to 90.08% for TRIC. The traits under a strong genotypic effect were FL, CE, and CI, exhibiting TSS values above 80%. On average, G accounted for 43.01%, 75.28%, and 60.94% for yield-related traits, fruit characteristics, and chemical features, respectively. The partitioning of the TSS% in the other components affecting the variation indicated that Y and G × Y accounted for, on average, 16.47% and 21.77%, respectively, resulting in, generally, less influence of Y and G × Y than G on the analyzed parameters.
Among the assessed traits, GY and RY were the traits most affected by seasonal fluctuations, showing a TSS of 65.21% and 37.72%, respectively, while for all fruit traits, except puffiness, the effect was very little, with a variation of 1.96% on average. The G × Y interaction was, on average, higher for yield traits (TSS equal to 25.91%) with respect to fruit characteristics (TSS equal to 20.99%) and their main chemical features (TSS equal to 18.46%). Among these, HI and PER were those for which the variation was highly influenced by the G × Y interaction, with values of 33.30% and 36.02%, respectively. Table 4. Analysis of variance and significant levels for the genotypic (G) and environmental effects due to growing season (Y) for Solanum habrochaites (SH) ILs and the SL control, and the combined effects for the traits evaluated.
The index of phenotypic plasticity (CRP) ( Table S4) was highly variable for the traits tested in the 39 ILs and the control, ranging from 0.05 (UM) to 5.40 (FW). On average, the fruit shape and color traits were more stable, with average plasticity values below 0.62. The yield traits, instead, exhibited higher plasticity, with average indexes above 0.79. In the genotypes tested, most of the values were lower than 1.00. Several ILs had CRP less than 1.00 for all the assessed traits, such as TA1128, TA1536, TA1537, and TA1312. We instead found values above 2.00 for TA1105 on chromosome 2 and for several lines on chromosomes 2, 4, and 5. Only TA1316 (chromosome 8) had CRP values < 0.93 for all traits. The yield traits exhibited high plasticity for TA1276 (HI = 4.68), TA1542 (PW = 4.16), and TA1562 (UM = 4.88), while the maximum value (5.40) was found for FW in TA1548. Instead, the fruit shape and color traits had lower values of plasticity, with an average index below 0.61.

Multivariate Analysis and Correlations between Traits
The principal component analysis (PCA) in the first two dimensions explained 50.94% of the total variance ( Figure 2).  Table 3.
The genotypes were evenly distributed on the two axes of the biplot. The first component accounted for 27.95% of the total variance and was positively correlated with all the traits except SSC and BY. The second component accounted for 22.98% of the variance, being positively correlated with all yield-related traits and negatively correlated with the remaining ones. For the entire set analyzed, the total phenotypic variation was explained by the first 19 components (data not shown). The yield-related traits explained the largest part of the variation in the first two components, showing a contribution of the variables with values above 12% (Table S5). The projection of genotypes on the two-dimensional PCA graph showed the wide range of phenotypic variations enclosed in the SH IL set, highlighting how wild alleles can improve diverse traits. For each growing season, correlograms between pairs of variables were generated using Bonferroni correction (p < 0.01) (Figure 3). Correlations mainly occurred among the same trait categories. In all cultivation trials, strong negative correlations were found for UM and GY, PW and HI, and PW and UM, while strong positive correlations were found between FD and FL, FW and FD, CE and CI, and RY and HY, as well as among BY, BM, TY, and RY. Only in the first and second years of cultivation, TRIC was correlated negatively to UM and positively to GY ( Figure  3a,b), while, during the third year, it was negatively correlated with BM (Figure 3c). Furthermore, a significant positive correlation between FW and PER was detected in the third The control is indicated with a square in red font. The direction and distance from the center of the biplot indicate how each trait contributes to the first two components. The explanation for trait acronyms is reported in Table 3.
The genotypes were evenly distributed on the two axes of the biplot. The first component accounted for 27.95% of the total variance and was positively correlated with all the traits except SSC and BY. The second component accounted for 22.98% of the variance, being positively correlated with all yield-related traits and negatively correlated with the remaining ones. For the entire set analyzed, the total phenotypic variation was explained by the first 19 components (data not shown). The yield-related traits explained the largest part of the variation in the first two components, showing a contribution of the variables with values above 12% (Table S5). The projection of genotypes on the two-dimensional PCA graph showed the wide range of phenotypic variations enclosed in the SH IL set, highlighting how wild alleles can improve diverse traits. For each growing season, correlograms between pairs of variables were generated using Bonferroni correction (p < 0.01) (Figure 3). Correlations mainly occurred among the same trait categories. In all cultivation trials, strong negative correlations were found for UM and GY, PW and HI, and PW and UM, while strong positive correlations were found between FD and FL, FW and FD, CE and CI, and RY and HY, as well as among BY, BM, TY, and RY. Only in the first and second years of cultivation, TRIC was correlated negatively to UM and positively to GY (Figure 3a,b), while, during the third year, it was negatively correlated with BM ( Figure 3c). Furthermore, a significant positive correlation between FW and PER was detected in the third season. In all trials, the soluble solids were negatively correlated with the yield traits and harvest index, although a significance was only observed during the last year. The degree of relationship between the elements of the phenotypic matrices revealed significant correlations between the data of the independent growing seasons, highlighting the strongest correlations between the second and third seasons of the evaluation: Y1/Y2 (r = 0.2, p < 0.0001), Y1/Y3 (r = 0.361, p < 0.0001), and Y2/Y3 (r = 0.607, p < 0.0001). season. In all trials, the soluble solids were negatively correlated with the yield traits and harvest index, although a significance was only observed during the last year. The degree of relationship between the elements of the phenotypic matrices revealed significant correlations between the data of the independent growing seasons, highlighting the strongest correlations between the second and third seasons of the evaluation: Y1/Y2 (r = 0.2, p < 0.0001), Y1/Y3 (r = 0.361, p < 0.0001), and Y2/Y3 (r = 0.607, p < 0.0001).  Table 3.

Genetic Regions Underlying the Phenotypic Variations
The overlaps between the introgressed fragments of the studied ILs determined the chromosomal regions (or BIN) carrying the wild alleles into the cultivated background. The size of each BIN was determined considering the overlapping regions, plus the halfintervals between the flanking markers and nearest mapped loci. The markers' positions were updated based on the Tomato-EXPEN reference map [32]. In total, 61 BINs were detected, ranging from a minimum of three on chromosomes 2, 5, 7, 8, 9, and 10 to a maximum of seven on chromosome 1 (Figure 4). Based on the mean comparisons of each ILs with the recurrent parent (Table S2), we detected 75 associations for all traits except for puffiness and green yield. Forty-two associated regions were found across two independent seasons, whereas 33 were across all trials, covering 17 genomic regions on nine chromosomes and including sixteen out of the nineteen evaluated traits. Except for chromosome 8, wild alleles were responsible for the quantitative variations for the studied traits in two or three trials (Table 5).

Figure 3.
Pearson's rank correlation coefficients between the pairs of phenotypes. The correlation coefficients are indicated in each cell. Colored correlations are those with p-values < 0.01 after Bonferroni correction. The color intensity is directly proportional to the coefficients. On the right side of the correlogram, the color legend shows the correlation coefficients and the corresponding colors. (a) The correlogram for the traits scored in year 1, (b) year 2, and (c) year 3. The trait acronyms are in Table 3.

Genetic Regions Underlying the Phenotypic Variations
The overlaps between the introgressed fragments of the studied ILs determined the chromosomal regions (or BIN) carrying the wild alleles into the cultivated background. The size of each BIN was determined considering the overlapping regions, plus the halfintervals between the flanking markers and nearest mapped loci. The markers' positions were updated based on the Tomato-EXPEN reference map [32]. In total, 61 BINs were detected, ranging from a minimum of three on chromosomes 2, 5, 7, 8, 9, and 10 to a maximum of seven on chromosome 1 (Figure 4). Based on the mean comparisons of each ILs with the recurrent parent (Table S2), we detected 75 associations for all traits except for puffiness and green yield. Forty-two associated regions were found across two independent seasons, whereas 33 were across all trials, covering 17 genomic regions on nine chromosomes and including sixteen out of the nineteen evaluated traits. Except for chromosome 8, wild alleles were responsible for the quantitative variations for the studied traits in two or three trials (Table 5).   Table 5. List of putative chromosome regions underlying the variations of the traits (QTL) identified across two or three independent trials. PV = percent of the phenotypic variance estimated and -= not available. In bold are indicated the positive PV (%). In capital letters are indicated those QTLs that increased the performance of the traits compared to the control. BIN sizes are determined considering the overlapping regions plus the half-intervals between flanking markers and the nearest mapped loci (see Figure 4).

Trait
ILs

Yield Traits
For the total yield and red yield, in all cases, the wild alleles were responsible for the reduction of productivity, with the PV ranging from −52.35% to −33.51% for TY and from −52.30% to −41.28% for RY. In total, six QTLs were found to be significant in two growing seasons, while three robust associations were found across the three years, being located in both the short and long arm of chromosome 1 (ty1.2 and ry1.3) and on chromosome 10 (ry10.1).
For the uniformity of maturity, two regions on chromosomes 1 and 4 were responsible for a decrease of the traits, with a PV −15% in both instances. The wild alleles exerted the highest effect on chromosome 1, reducing the trait across all the seasons.
As concerns the biomass and harvest index, a total of nine QTLs were detected, with the larger effect of the wild alleles on the former trait leading to a reduction ranging from −52.03% to −40.24%, whereas, for the latter, the PV ranged from −25.52% to −12.14%. Two QTLs had significant effects during the three consecutive seasons, being responsible for the reduction of the traits for −40.24% (bm1.2) and −25.52% (hi3.4). On the contrary, only for the plant weights, the wild alleles located in a 15-cM chromosomal region on chromosome 3 were responsible for an increase (PW3.6), with a PV of 73.92%.

Fruit Traits
Numerous alleles were found to be responsible for the variations of the fruit traits. The fruit weight was affected by eight significant QTLs; in all cases, the SH alleles reduced the fruit sizes. The genomic regions underlying the FW were identified on seven chromosomes. Five QTLs were robust across the three years, with the strongest effect found for fw11.1, which accounted for a percentage of the phenotypic variance lower than −43.14%. The lowest effect was found for fw3.5, with a PV of −27.91%. Three QTLs were identified in two seasons on chromosomes 9, 10, and 11, with PV% ranging from −35.63% (fw10.3) to −32.01% (fw11.3).
For fruit morphological traits, a total of nine and seven QTLs were found to be significantly associated with the fruit lengths and fruit widths, respectively. The SH alleles decreased the two dimensions in all cases, leading to fruits with small sizes. Robust QTLs were found across the three years. A robust cluster was found on chromosome 7 (fl7.3 and fd7.3), with a PV of −19% in both instances. Closely located QTLs were found on chromosome 3 for both traits (fl3.4 and fd3.5) and on the bottom part of chromosome 12 for the fruit widths (fl12.3 and fd12.4). Overall, the reduction of the lengths and the widths of the fruits were lower than 20% for all the detected genomic regions.
In total, five QTLs for the fruit shape were detected; for three identified across all the trials, the SH alleles caused the fruit to be more flattened (decreasing of the FS). Among these, fs3.1 and fs11.1 had the highest and the lowest effects, respectively, with PV% ranging from −12.73% to −5.96%. On the bottom of chromosome 1, we found alleles leading to more elongated fruits, with increased shapes of 7.42% (FS1.5) and 7.28% (FS1.4) during seasons two and three, respectively.
The wild SH alleles were responsible for both decreasing and increasing the intensity of the fruit colors. Four QTLs for the internal and external colors of the fruits were found to colocalize at the bottom of chromosomes 1 and 2. In all instances, the wild alleles led to a reduction of color intensity, from red to orange/yellow. The strongest effects were observed on chromosome 2, with a PV of −53.95% for the internal color across the three years of evaluation. Minor effects were instead found on chromosome 1, with variations on growing seasons two and three for the internal (ic1.4) and external (ec1.4) colors, respectively. The QTLs involved in the increase of the fruit color intensity were detected across two seasons on the long arms of chromosomes 4, 9, and 11, with a PV ranging from 18.97 (EC11.3) to 25.05 (IC4.5). Although not colocalizing, neighboring QTLs were found at the base of chromosome 4 in a region spanning 96.75-119.0 cM.
Finally, the QTLs involved in pericarp thickness were identified on chromosomes 5 and 11, respectively. The former was detected across all seasons, increasing the pericarp thickness up to 32.49%, while the latter was detected only in two trials, leading to a reduction of −38.69%.

Chemical Traits
For the chemical traits, we found 11 associated regions in two or three growing seasons. In all cases, the SH alleles increased the soluble solids with a PV in the range of 24-27%. Most of the significant QTLs were detected during two seasons. Only SSC1.5, on the long arm of chromosome 1, had the greatest overall effect, being detected across all seasons. Despite the increase of the soluble solids, the product brix × red yield decreased in percentage, ranging from −57.73% to −39.33%. In tomatoes, this trait provides the estimation of the paste amount during processing. We also found four QTLs associated to a reduction of pH on the long arms of chromosomes 6 and 9 in two seasons (average PṼ 4.5%) and on chromosomes 1 and 11 across three seasons, with an average PV of 8%.

Trichomes
The amount of trichomes is of interest, as it provides an estimate of the level of tolerance/resistance to virus-carrying insects such as aphids and thrips [42]. This protection mechanism is exercised both as a mechanical opposition and through the production of volatile repellent substances contained in the trichromatic glands [26]. Although there are very advanced methods of investigation based on microscopic techniques, the scale we used serves to give an estimate of the putative QTLs underlying the trait. Three highly significant QTLs were found on chromosomes 4, 7, and 10 and were related to TA1280, TA1304, and TA1551, respectively. In all cases, the wild alleles provided an increase of the trichomes number in comparison to the control in a percentage range of 100-200%.

Performance of S. habrochaites Introgression Lines
Changes of the agricultural scenario in the coming decades will require efforts from breeders in finding novel sources of variations for the improvement of crops [7]. Wild tomato species have evolved in a wide range of habitats and, as a result of natural selection, represent a useful reservoir of genes for various agro-morphological, physio-biochemical, and resistances traits [2]. This exotic genome, once introgressed in the cultivated background, can be successfully exploited, breaking up the association with undesirable traits (linkage drag) and overcoming any occurring breeding barriers [43]. In the present study, a comprehensive approach aimed at assessing agronomic qualitative performances was carried out in multi-season trials on a set of 39 S. habrochaites LA1777 ILs. We firstly examined the performance of the ILs by focusing on the genetic and environmental factors affecting trait variations. Subsequently, we investigated those genomic regions underlying the studied traits and identified putative QTLs. A wide range of diversity was found among ILs for the considered traits with a medium-high level of heritability, which was maximized for fruit shape features. Despite the presence of genetically diverse lines and of additional introgressions for 26% of the lines, the fruit traits were minimally affected by environmental factors. In the set analyzed, we found trait plasticity depending on the chromosomal regions and type of traits, due likely to the different types of gene regions interacting with the environment.
Overall, the whole set of introgression lines resulted in being very promising for soluble solid contents rather than yield traits. The level of soluble solids measured as brix degrees is a measure of the total content of sugars and, in a minor extent, of the organic acids, together with small amounts of dissolved vitamins, fructans, proteins, pigments, phenolics, and minerals [44]. All these compounds influence the flavor and sweetness, which are important precursors of sensory quality [45].
Four lines at the bottom of chromosomes 1 (TA523 and TA1223) and 4 (TA1473 and TA1475) were previously used as a basis to fine map the agronomic traits [10,11]. Our results agree with earlier observations reporting an increase of soluble solids and a decrease of yield traits at the bottom of chromosome 1 in the overlapping region between TA523 and TA1223 [10]. Moreover, we confirmed a general decrease in fruit weight and an increase in the soluble solids content and color intensity at the bottom of chromosome 4 [11]. Regarding productivity, the differences in performances of the control genotype and the IL set in each cultivation year were almost stable. We observed a decrease of the yield and uniformity of maturity from the first to the last season. As expected, similar trends of correlations were found across the years, particularly within the same category of traits. Interestingly, significative negative correlations were found between the yield traits and soluble solids in the last growing season characterized by less rainfall, humidity, and average temperatures. It is likely that the combination of these effects led to a general decrease of TY and RY compared to the previous years in both E6203 and the ILs. Instead, the first season was characterized by the highest average temperatures and humidity, leading to the best yield values.
The soluble solids tended to be stable across the three seasons, although slightly higher in the less productive ones. Within the ILs, the contrasting values observed in the second year mainly concerned lines with multiple introgressions and/or low brix values in the other seasons, suggesting a major effect of the environment in these cases. Previous evidence in tomatoes reported how increasing the water stress led to the increase of the soluble solids content [46]. Moreover, the rainfall and temperature changes could significantly impact the overall performance of the tomatoes [47].
Although the used drip irrigation restored the entire water requirements at the root level, it is likely that the whole soil moisture was affected by the conditions that occurred in the last season, leading to the observed correlations. Furthermore, it is known that S. habrochaites is a cold-tolerant species; therefore, it is possible to suggest that lower temperatures may have a role in the increase of the soluble solids. Since no previous evidence was reported in S. habrochaites acc. LA1777, further investigation is required.

QTL Mapping
Despite S. habrochaites being recognized as an important source of alleles of agronomic interest, the potential of the existing exotic libraries has not yet been extensively exploited. We identified 75 regions likely with underlying trait variations, and the multi-evaluation trials allowed the identification of 33 robust QTLs, removing any possibility of false discoveries.
Our results agree with previous studies involving diverse types of exotic libraries and experimental mapping populations developed from wild relatives [8,29,30,[48][49][50]. Many QTLs reducing the yield performances and fruit weights confirm the preponderant effect of the wild Solanum alleles on these traits. We found specific regions associated with an increase in plant weights, soluble solids contents, internal and external color intensities, pericarp thickness, and trichome density. The results of this study confirmed that S. habrochaites could be a potential source of favorable alleles for horticultural traits, despite the inferior phenotype compared to the cultivated parent. Some genomic regions seem to be conserved across wild tomato relatives; as an example, the increase of soluble solids at the bottom of chromosome 1 fall in a region previously reported carrying alleles improving the trait in other wild species, including S. pennellii, S. neorickii, and S. peruvianum [8,[48][49][50]. Furthermore, the additional QTLs detected on chromosomes 3 (SSC3.6) and 10 (SSC10.1) were also reported in S. pennellii introgression lines [8] and, for the one on chromosome 10, in advanced backcross lines of S. peruvianum [48,50]. These results confirm the existence of shared regions in various tomato wild species carrying alleles for improving the quality of cultivated ones. Clusters of QTLs involved in fruit weight and morphology were found at the long arms of chromosomes 3, 7, 11, and 12, suggesting a gene linkage. For all, the underlying alleles led to a reduction of fruit weights and shapes in agreement with the wild phenotype, which has an average weight fifty-fold less than the recurrent parent and rounder fruits. Fruit color is one of the main attributes of quality and an indicator of the contents of health-related compounds. Diverse orange, yellow, pink, and black tomatoes are becoming increasingly popular by consumers who appreciate a rich and varied diet. The decrease of internal and external fruit color at the bottom of chromosome 2 was associated with fw2.2, a major gene modulating fruit sizes [33]. The same gene was reported being associated with increased color in S. pimpinellifolium due to a pleiotropic effect between the decreasing fruit size and pigment concentration [51]. Pleiotropy between the fruit size and color occurs also in SH, although with the opposite effect. In S. pennellii, a decrease in color was shown by the mutation of yellow flesh (r) on chromosome 3, responsible for the loss of function of the phytoene synthase, a key enzyme of the carotenoid biosynthetic pathway [52]. The mutation was carried by IL3-2, which defines the same region of TA1539 and TA1541. However, these two lines had additional introgressions, which likely masked a possible effect of any occurring mutation.
Therefore, the additional alleles part of the SH genomic background could be involved in fruit pigmentation. A strong QTL responsible for the increase of trichome density was found on chromosome 10. Recently, Barrantes et al. [15] identified a major QTL in the same chromosomal region in S. pimpinellifolium ILs. In agreement with our study, the trait was scored by the same scale. Trichomes are associated with insect resistance, although no candidate gene has been reported on chromosome 10. These results suggest that further studies are needed for a better investigation of the underlying region, although it cannot be excluded that this increase in the density of trichomes does not correspond to a resistance.

Future Prospects for Using S. habrochaites ILs in Tomato Breeding
Beyond QTL identification, the advantage of using introgression lines relies on the possibility of their direct use in breeding programs as parent lines. This offers an invaluable opportunity for transferring new traits for the improvement of modern varieties. The set of SH ILs analyzed in this study resulted in being promising for the diverse phenotyped traits. Furthermore, the number of QTLs detected in specific lines facilitates gene pyramiding; as an example, the improvement of soluble solids and pericarp thickness could be obtained by crossing TA1287 either with TA1223 or TA523, which, in addition, showed a slight increase of pH. The synergy of these QTLs could improve the processing yields, with less acidic fruits with greater pulp consistency. For a successful introgression program, it is very important to choose phenotypes stable across a wide range of environmental conditions. To that end, QTLs with large genetic effects are preferred. However, it must be considered, the occurrence of constraints such as pleiotropy, linkage-drag, epistatic effects, and relationships with other beneficial traits. For example, introgressing soluble solids is often challenging due to inverse correlations with the yield [53]. Moreover, pleiotropic effects may occur between the fruit size and color [51].
Despite one of the advantages in holding only a fraction of the entire wild genome being the reduction of fertility problems [8], it must be recognized that several S. habrochaites ILs still carry multiple wild segments, as well as large introgressed fragments, covering, in some instances, more than half the chromosome. Therefore, the linkage between the beneficial genes introduced during the backcrossing with the deleterious ones (linkage drag) could reduce the potentiality of ILs for both breeding and QTL mapping purposes. In the current set, we demonstrated how the availability of lines with reduced introgressions led to a major accuracy in the identification and positioning of QTLs with positive effects (e.g., SSC1.5 and PER5.1). While linkage drag occurs for large introgression lines such as TA1304, with both the increasing of the trichome density and reduction of the fruit weight and size.
Our results suggest, in some instances, the necessity of generating new sub-ILs able to break up the linkage occurring with undesired traits and facilitate the fine mapping of QTLs. To that end, genomic platforms such as SNP (single nucleotide polymorphism) arrays or reduced-representation sequencing (e.g., ddRAD-seq: Double digest restrictionsite associated DNA and GBS: Genotyping by sequencing) could serve as a key step for the generation of high-density marker data, speeding up the deep exploitation of the current S. habrochaites exotic library. These tools have been successfully applied in tomatoes and other related species, providing large-scale and cost-effective genome scans and generating novel knowledge for genomic-assisted breeding [54,55].

Conclusions
The present study aims to provide a broad overview of S. habrochaites introgression lines in relation to the morpho-agronomic performance and chemical properties. Through the investigation of phenotypic performances, the G × Y interaction, and the QTL analysis, we provided novel information for genetic and breeding purposes. This study represents a step toward the exploitation of the hidden properties of S. habrochaites, to be implemented through the development of near-isogenic lines assisted by genomics approaches.
Supplementary Materials: The following are available online at https://www.mdpi.com/2073-4 395/11/1/38/s1: Table S1: Parameters describing the trait performances of the recurrent parent "E6203" and the 39 SH ILs for each growing season. For each trait, the significant differences between the ILs and the recurrent parent is reported. Means in the same column with different letters are significantly different, according to Tukey's honest significant difference test; Table S2: Mean values for the introgression lines and the control genotype in each growing season (Y1, Y2, and Y3). In bold and with asterisks are indicated those means that significantly differ compared with the control, according to Dunnett's test; Table S3: Percent of phenotypic variation with respect to the control for the introgression lines in each growing season (Y1, Y2, and Y3). In bold are indicated those values corresponding to significant differences according to Dunnett's test; Table S4: The phenotypic plasticity index (PPI) for the traits analyzed in 39 S. habrochaites ILs and the control S. lycopersicum E6203. PPI values > 1.0 are in bold and gray highlighted. For each IL is indicated the chromosomes holding introgressions; Table S5: Variable contribution (VarPC), correlation coefficient (CorrPC), and Eigenvectors (Egv) for fruit descriptors in the two first principal components. Highlighted traits explained a variation >12% in the first or second PC.