An F₁ population of 91 individuals derived from the cross between ‘Marcona’ and ‘Marinada’ (MC × MI) was used in this study. ‘Marcona’ is a highly valued traditional Spanish cultivar, known for its characteristic rounded kernel shape and high level of fatty acids (Martín Carratalá et al., 1998; Calle et al., 2024). ‘Marinada’ (‘Lauranne’ × ‘Glorieta’) is a self-compatible breeding cultivar released by IRTA in 2008, noted for its very sweet kernel and high percentage of soluble sugars (Vargas et al., 2008; Romero et al., 2011). Seedlings were grafted onto ‘Garnem^®’ rootstock and planted at 4 × 1.8 m in 2015. The MC × MI population and its parents are maintained at the IRTA Mas Bové experimental station (41.170723 N, 1.172942 E) under standard agricultural practices.

The MC × MI population and both parents were evaluated for several physical and chemical traits. The physical traits included kernel weight, crack-out percentage, kernel size (length, width, and thickness), kernel shape (roundness, globosity, shoulder, and symmetry), and tegument and kernel color (L*, a*, b*). The kernel chemical traits assessed were kernel protein, fiber and fat content, and fatty acids profile (myristic, palmitic, palmitoleic, margaric, cis-10-heptadecenoic, stearic, oleic, vaccenic, linoelaidic, linoleic, arachidic, cis-11-eicosenoic, and cis-11,14-eicosadienoic).

From 2018 to 2021, fifty mature fruits were randomly collected from each individual of the F₁ MC × MI population and its parents. Fruits were considered mature when the mesocarp was fully dry, split along the fruit suture, and the peduncle was near complete abscission. After dehulling, nut weight was measured, shells were cracked, and kernel weight (KWe) was recorded using an electronic scale. Crack-out percentage, calculated between 2019 and 2021, was determined as the ratio of kernel and nut weights. After that, kernel length (KLen), width (KWidth), and thickness (KThick) was measured with a digital Vernier caliper. Kernel roundness (KRound) and globosity (KGlob) were estimated using the ratios width/length and width/thickness, respectively. Kernel size and shape traits were also assessed using image analysis. For that, a standard photo of six kernels was taken per individual ( Figure 1 ), and images were analyzed using the Shape Analyzer (SA) software (Jurado, 2024). This deep learning-based tool automatically detects almond kernels from images and measures them. The parameters obtained included kernel length (KLen_SA), width (KWidth_SA), roundness (KRound_SA), symmetry based on ‘structural similarity index measure’ (SSIM) (Sym_SSIM), and Jaccard index (Sym_Jacc). Both indices determine kernel symmetry by comparing the similarity between the two sides of the kernel along the longitudinal axes using digital images ( Figure 1B ). For that, the box containing the kernel ( Figure 1B ) was divided in half and right-flipped to compare the similarity between both sides using the SSIM and Jaccard coefficient metrics, with a value of 1 indicating perfect symmetry. These traits were phenotyped in 2020 and 2021. Additionally, in 2020, the kernel shoulder, referring to anomalies in the basal area of the kernels ( Figure 1 ), was visually assessed using a scale from one (no shoulder) to five (marked basal anomaly as observed in Figure 1C ; KShoul_Vis), and measured using the Tomato analyzer software (Gonzalo et al., 2009) to determine the angle of the kernel shoulder (KShoul_Angle) as shown in Figure 1 . Tegument (the kernel skin) and kernel color were determined with a Minolta Chroma Meter tri-stimulus color analyzer (CR-3500D; Minolta, Ramsey, NJ, USA) calibrated to a white porcelain reference plate using a CIELAB scale with color space coordinates L*, a* and b* (Tcolor_L, Tcolor_a, Tcolor_b, Kcolor_L, Kcolor_a, Kcolor_b). Tegument and kernel color were measured from 2018 to 2021.

Fat content was analyzed by Soxhlet method, using 5 – 6 g of ground blanched almonds and petroleum ether (boiling point 40 to 60 °C) for 7 h in Soxhlet apparatus (Romero et al., 2021; Etheridge et al., 1998). The fat content was expressed as a percentage (%), calculated based on the weight of the extracted fat relative to the total weight of the sample. This trait was measured only in 2021.

Crude protein was analyzed from 2018 to 2021 by Dumas combustion procedure using Leco FP-528 analyzer as described in Romero et al. (2021). Briefly, 0.2 g of grounded sample was weighed in a porcelain sample holder (boat) for introduction into the combustion chamber (850 ± 1 °C) utilizing an automated sample loader. The combustion process converts covalently bound 130 nitrogen into nitrogen gas (N₂) that is quantified by passing the gas through a conductivity cell and converted to protein by multiplying a factor of 6.25. The results were expressed as a percentage (%) of protein content.

Crude fiber was measured from 2019 to 2021 using 1 g of the grounded sample. The sample was treated with boiling 0.26 N sulfuric acid for 30 min, followed by boiling 0.23 N potassium hydroxide for another 30 min (Romero et al., 2021). The extracted residue was dried at 103 ± 1 °C for 3 h, weighed, then placed in a furnace (550 ± 1 °C for 3 h), and finally, the ashes were weighed. These results were expressed as a percentage (%) of ash content.

Fatty acids were analyzed by gas-chromatography with a flame ionization detector (GC-FID) using a capillary column (Romero et al., 2021). The fatty acid methyl esters (FAMEs) were prepared by trans-esterification with 0.5 M potassium hydroxide, following the official method UNE-EN ISO 5509:2000. FAMEs (1 mL) were separated using a gas-chromatograph (HP 6890; Agilent Technologies, Barcelona, Spain) equipped with an FID detector and a capillary column [30 m · 0.25 mm i.d. (HP-Innowax, Agilent Technologies)]. The carrier gas was helium, with a flow rate of 1 mL/min. The injector and detector temperatures were 220 and 275°C, respectively. The FAME identification was based on retention time relative to those of a standard FAME mixture (Sigma-Aldrich, Madrid, Spain). Fatty acids were measured only in 2021. Results were expressed in mg/g as average values of three replicates.

In this study, LSmean, which provides a way to obtain mean trait values that are adjusted for fixed effects, was considered to predict the effect of the year, assuming that the environmental conditions of every year affect all individuals in the population in the same way. Therefore, for traits with more than one year of data, the LSmean was calculated according to the following equation:

where P i j is the phenotypic value of the i-th individual in the j-th year, β 0 is the intercept (overall mean phenotypic value), β i   is the effect of the i-th individual, β y e a r j is the effect of the j-th year, and e i j is the error term for the i-th individual in the j-th year. The coefficient of determination (R²) was used to measure how the LSmean fits the data and how well it can predict outcomes. R² was determined using the following formula:

where SST is the sum of squares of the residual errors, and SST is the total sum of the errors. R² ranges between 0 and 1, with values closer to 1 indicating better LSmean prediction.

Pairwise correlation coefficients were calculated for all traits (single-year data and LSmean data) using the JMP^® software (JMP^®, Version 16, SAS Institute Inc.) Significance was calculated using the Spearman correlation coefficient significant (p<0.001). Normal distribution was assessed using the Shapiro-Wilk test (p<0.05) using JMP^®.

Total genomic DNA from the 91 individuals and two parents was isolated from young leaves using the CTAB method (Doyle and Doyle, 1990), adapted for 96-well plates. The DNA samples were genotyped with 33 SSRs and the almond Axiom™ 60K SNP array (Duval et al., 2023).

For the SSRs, a set of markers that were heterozygous in ‘Marcona’ or in ‘Marinada’ was selected. In genomic regions not well covered with the initial set of SSRs, new SSRs were designed using the almond reference genome (Alioto et al., 2020). Primers were designed using Primer 3 (http://primer3.ut.ee, v4.1.0) with default parameters. The list of SSR markers used is presented in Supplementary Table 1 . PCR reactions were conducted in a final volume of 10 μL containing 200 ng of genomic DNA, 1x NH₄ reaction buffer, 1.5 mM MgCl₂, 0.2 mM dNTPs (10mM), 0.2 μM of each marker, 1 U of BIOTaq (Bioline, London, UK) and HPLC H₂O to reach the final volume. PCRs were performed in a GeneAmp PCR System 9700 thermal cycler (Applied Biosystems, CA, USA) with the following conditions: initial denaturation at 94°C for 1 min, 35 cycles of denaturation at 94°C for 15 s, primer annealing at a specific temperature for each primer for 15 s, extension at 72°C for 30 s, and a final extension at 72°C for 5 min. Forward primers were designed with a generic fluorochrome sequence at the 5’ ends (FAM, VIC, NED, or PET). PCR products were added to 12 μL of deionized formamide containing 0.35 μL of GeneScan500 LIZ size standard (Applied Biosystems, CA, USA) and heated at 94°C for 3 min. Capillary electrophoresis was performed in an ABI Prism 3130xl automated sequencer (Applied Biosystems, CA, USA). GeneMapper v5.0 software (Applied Biosystems) was used for SSR allele sizing.

The almond Axiom™ 60K SNP array genotyping was performed on an Axiom GeneTitan™ (ThermoFisher Scientific) system at the INRAE Gentyane platform in Clermont-Ferrand (France). Genotypic data was retrieved using the Axiom Analysis Suite (ThermoFisher, 2017). Samples were filtered following the Axiom best practices workflow, setting an average call rate > 95. We then filtered out SNPs with the following characteristics: (i) monomorphic SNPs in the progeny; (ii) heterozygous SNPs in ‘Marcona’ and ‘Marinada’ but with only two genotypic classes in the progeny; (iii) homozygous SNPs in ‘Marcona’ and heterozygous in ‘Marinada’ with three genotypic classes in the progeny; (iv) homozygous SNPs in ‘Marinada’ and heterozygous in ‘Marcona’, but with three genotypic classes in the progeny. After filtering based on segregation, missing data, and putative genotyping error identification and imputation were performed with AlphaFamImputed software (Whalen et al., 2020) using default settings. We then ordered the SNPs based on their physical position in the ‘Texas’ almond genome v2.0 (Alioto et al., 2020) and phased them manually. Finally, a set of bins (i.e. groups of SNPs with identical genotypes for all the individuals) was established, with each bin separated from the adjacent bin by at least one recombination event.

Finally, three linkage maps were built using JoinMap 5^® (van Ooijen, 2018). The MC × MI map was constructed using all the previously selected bins, the ‘Marcona’ map was built using only bins heterozygous in ‘Marcona’, and the ‘Marinada’ map was built using only bins heterozygous in ‘Marinada’. For map construction, a minimum logarithm of odds (LOD) score of 10 was selected for SNP grouping. Makers showing segregation distortion higher than 0.01 were excluded from linkage mapping unless other distorted markers with similar ratios surrounded them. The maximum likelihood algorithm with default parameters and the Kosambi mapping function were used (Kosambi, 1944).

LSmean and single-year data were analyzed for QTL detection in the three linkage maps (‘Marcona’, ‘Marinada’, and MC × MI) using MapQTL 6.0^® (van Ooijen, 2009). Interval mapping (Lander and Botstein, 1989) and multiple QTL mapping (Jansen and Stam, 1994) strategies were used for QTL discovery. To determine the significance threshold for each QTL, LOD was calculated for each linkage group (LG) and trait using the permutation test (1000 permutations) at a 95% significance level (p < 0.05) (Lander and Botstein, 1989; van Ooijen, 1992). As most traits had a significance level similar to 3, a final threshold of LOD 3 was used in all cases. The QTL confidence interval was defined as the LOD - 1. For traits with non-normally distributed data, we used the Kruskal-Wallis non-parametric test in MapQTL, retaining only QTLs that were significant in this test (p< 0.01) and also significant with interval mapping (LOD > 3.0). QTLs were named according to the recommendations for standard QTL nomenclature and reporting of the Genome Database for Rosaceae (Jung et al., 2019). Graphical representation of linkage maps and QTLs was created using MapChart (Voorrips, 2002). The physical positions of the QTLs are based on the almond reference genome of ‘Texas’ v3.0-F1 (Castanera et al., 2024).

We obtained phenotypic data for 22 traits, including both kernel physical and chemical characteristics. The physical traits measured were kernel weight, crack-out percentage (crack-out), kernel size (length, width, and thickness), kernel shape (roundness, globosity, shoulder, and symmetry), and tegument and kernel color (L*, a*, b*). The kernel chemical traits included kernel protein, fiber, and fat content. For the fatty acid profile, we measured myristic, palmitic, palmitoleic, margaric, cis-10-heptadecenoic, stearic, oleic, vaccenic, linoelaidic, linoleic, arachidic, cis-11-eicosenoic, and cis-11,14-eicosadienoic acids.

A total of 9243 SNPs, representing 15.3% of the almond Axiom™ 60K SNP array, were obtained from the SNP filtering procedures and used for linkage mapping along with 33 simple-sequence repeat (SSR) markers. Of these SNPs, 3810 were heterozygous in ‘Marcona’ and homozygous in ‘Marinada’, 4060 were heterozygous in ‘Marinada’ and homozygous in ‘Marcona’, and the remaining 1373 SNPs were heterozygous in both parents. Additionally, 33 SSR markers were included.

For all constructed linkage maps (‘Marcona’, ‘Marinada’, and ‘MC × MI’), markers were grouped into eight LGs, each corresponding to an almond chromosome ( Supplementary Table 5 ). The ‘Marcona’ map covered 448.0 cM and presented 370 bins (i.e., groups of markers in unique genetic positions) ( Table 1 ). Similarly, the ‘Marinada’ map had 316 bins distributed along 447.0 cM. The ‘MC × MI’ showed higher marker saturation than parental maps, with 685 bins covering 540.8 cM ( Table 1 ). In the three linkage maps, LG1 covered the largest genetic distance and included the highest number of markers. The average distance between consecutive bins was less than 1 cM in the ‘MC × MI’ map, whereas for the parental maps, this distance ranged from 1.15 (‘Marcona’; LG1) to 1.61 cM (‘Marinada’; LG7) ( Table 1 ). Additionally, the linkage maps showed good genome coverage without major gaps: the average physical distance between consecutive markers was 196, 555, and 612 Kbp in the ‘MC × MI’, ‘Marcona’, and ‘Marinada’ maps, respectively.

QTL analyses were performed using single-year data for all traits, and for those traits with more than one year of data, LSmean values were also used. All QTLs are presented in Supplementary Table 6 . In most cases, LSmean data identified the most consistent QTLs, which were detected in at least two different years through single-year QTL mapping. To simplify data presentation, we present QTLs using LSmean and those for which only one year of phenotypic data was available ( Table 2 ).

In total, 51, 18, and 24 QTLs were detected using the MC × MI, ‘Marcona’, and ‘Marinada’ genetic maps, respectively. After excluding QTLs that cosegregate across different maps, a total of 53 unique QTLs were identified ( Table 2 ). Most QTLs detected with the parental maps were also identified with the MC × MI map, except qP-KWe3.1 (only found on ‘Marcona’ map) and qP-KWidth4.1 (only found on ‘Marinada’ map). Conversely, 20 QTLs (qP-KWidth7.1, qP-KRound6.1, qP-KGlob2.1, qP-KGlob3.1, qP-KLen_SA6.1, qP-KWidth_SA1.1, qP-KWidth_SA2.1, qP-KWidth_SA4.1, qP-KWidth_SA7.1, qP-KRound_SA6.1, qP-Sym_SSIM8.1, qP-Sym_Jacc8.1, qP-Protein2.1, qP-Protein3.1, qP-Protein5.1, qP-Protein7.1, qP-Fiber1.1, qP-Fat3.1, qP-Fat5.1, qP-VaccenicAc1.1) were identified on MC × MI map but not on the parental maps. All consistent QTLs identified with the MC × MI map are presented in Figure 3 .

At least one QTL was identified for all traits except for the color parameters Tegument_L, Tegument_a, Tegument_b, Kernel_L, and linoelaidic acid. For kernel weight, two regions were detected on LGs 3 and 4 ( Table 2 ), with the major QTL (qP-KWe4.1) found at the beginning of LG4, explaining 29.0% of the phenotypic variation (PV) in the MC × MI map and 26.1% in the ‘Marinada’ map. This QTL showed an additive effect of the ‘Marinada’ allele (+0.11 g) increasing kernel size in the population ( Table 2 ). A less significant weight QTL, only detected in the ‘Marcona’ map, was also identified on LG3 (qP-We3.1; PV=15.4%). For crack-out percentage, a highly significant QTL (qP-Crack2.1) located on LG2 was detected in the three maps with a LOD of 17.2 and explaining 59.7% of the PV in the MC × MI map ( Table 2 ; Figure 3 ).

For kernel size and shape traits, a major QTL was detected at the beginning of LG1 (11.4-19.2 cM) for kernel length (qP-KLen1.1) using the ‘Marinada’ and the MC × MI maps, explaining up to 32.1% of the PV. Similarly, another major QTL for kernel width (qP-KWidth1.1) was identified in the central part of the same LG1 (45.3-57.4 cM) in the ‘Marcona’ (LOD=3.8; PV=18.3%) and MC × MI (LOD=4.9; PV=22.9%) maps. Overlapping with qP-KLen1.1, other consistent QTLs were identified for kernel roundness (qP-KRound1.1; LOD=4.1 to 7.8; PV=19.6 to 33.9%) and kernel length_SA (qP-KLen_SA1.1; LOD=6.3 to 6.8; PV=28.5 to 31.2%). In the same region of qP-KWidth1.1, we identified other QTLs for kernel width SA (qP-KWidth_SA1.1), kernel symmetry SSIM (qP-Sym_SSIM1.1), kernel symmetry Jaccard (qP-Sym_Jacc1.1), and kernel shoulder (qP-KShoul_Vis1.1) ( Table 2 ; Figure 3 ). For kernel thickness and kernel shoulder angle, QTLs (qP-KThick1.1 and qP-KShoul_Angle1.1) were identified in positions near qP-KLen1.1; however, the LOD plots do not clearly reveal whether these represent the same QTL or two distinct QTLs located in very close proximity ( Figure 3 ). Other stable QTLs for kernel width (LOD = 3.3 to 3.9; PV = 15.8 to 18.7%) were also identified on LGs 4, 7 and 8 (qP-KWidth4.1, qP-KWidth7.1, and qP-KWidth8.1). Overlapping with qP-KWidth4.1, a consistent QTL for kernel thickness (qP-KThick4.1) was also identified in ‘Marinada’ (LOD=4.2; PV=20.0%) and MC × MI (LOD=5.3; PV=24.6%) maps ( Table 2 ; Figure 3 ). Additional regions with kernel globosity QTLs were reported on LGs 2, 3, and 4, each showing similar LOD scores (3.8-5.0) and PV percentages (18.1-23.3%). Comparison of detected QTLs for kernel length, width, and roundness using manual determinations and Shape Analyzer (SA) software measurements revealed QTLs in similar positions, with only a few QTLs detected exclusively by one method or the other (qP-KRound1.1, qP-Round_SA2.1, qP-Width_SA2.1, qP-KWidth_SA4.1, qP-KLen_SA6.1).

Regarding the kernel chemical composition, four QTLs for protein content were mapped on LGs 2, 3, 5, and 7 on the MC × MI map ( Table 2 ; Figure 3 ). Among them, qP-Protein2.1 showed the highest LOD scores (4.5) and explained the largest proportion of PV (21.2%). The other three additional protein content QTLs (qP-Protein3.1, qP-Protein5.1, and qP-Protein7.1) were detected with lower significance, having LOD scores below 4 and explaining no more than 19% of the phenotypic variance ( Table 2 ). For fat content, two QTLs were found on LGs 3 and 5 of the MC × MI map, both with the same significance (LOD=3.9) and explaining similar PV (19%) ( Table 2 ). Notably, the protein and fat content QTLs on LG5 (qP-Protein5.1 and qP-Fat5.1) overlapped in the same region ( Figure 3 ). For fiber content, a unique QTL mapped on LG1 (14.8-23.3 cM) of the MC × MI map was detected, associated with a LOD score of 3.4 and explained 16.2% of PV.

For the fatty acid content, QTLs were detected for all compounds except linoelaidic acid, on MC × MI map within the same region at the bottom of LG1 (70-85 cM). The QTL with the highest LOD score and PV explained for these acids was reported for margaric acid (LOD=8.2; PV=37.9%), whereas the mystiric acid QTL presented the lowest significance (LOD=3.8; PV=19.7%). On the ‘Marinada’ map, all the same QTLs for fatty acids reported in the MC × MI map were mapped to the same position of LG1. The ‘Marinada’ allele for these QTLs was associated with reductions in the mean values of all acids within the population ( Table 2 ). On the ‘Marcona’ map, only QTLs for palmitic acid, margaric acid and oleic acid were detected, all found on the same LG1 interval as for MC × MI and ‘Marinada’ maps, with PV not exceeding the 20% ( Table 2 ). For these QTLs, the ‘Marcona’ allele was associated with increasing in fatty acid content. This effect was notable on oleic acid, the main fatty acid in almonds, where the ‘Marcona’ allele resulted in an increase of 22 mg per gram of sample.

Finally, two QTLs (qP-KColor_a1.1 and qP-KColor_b2.1) associated with kernel color were mapped on the MC × MI and ‘Marinada’ maps ( Table 2 ; Figure 3 ). qP-KColor_a1.1, mapped on LG1, explained 25.3% of the PV in MC × MI map, whereas qP-KColor_b2.1, found on LG2 was associated with PV of 27.1% in MC × MI ( Table 2 ).

Most traits phenotyped over multiple years showed significant correlations across those years. This suggests that the phenotypic data obtained is of high quality and that most traits exhibit high heritability. These findings are consistent with other studies on traits like kernel shape (Martínez-Garcia et al., 2019), which facilitates inheritance studies and supports the implementation of marker-based selection approaches.

Identifying correlations between traits can provide valuable insights for more efficient breeding and selection. We found several interesting correlations between physical and chemical traits: for example, kernel weight is correlated with kernel length and width, but these two are not correlated with each other. Additionally, kernel roundness is correlated with kernel length and width, but not with kernel weight. Kernel thickness is correlated with kernel width and shape but not with kernel length, which opens the possibility of discovering kernel shapes different than those of the parents used. Our results indicate that kernel weight is more related to the size than the shape of the kernels, as reported in recent studies (Lipan et al., 2022). High correlations between morphological traits were also found previously in progenies derived from ‘Blanquerna’ × ‘Vivot’ (Fernández i Martí et al., 2013) and ‘Nonpareil’ × ‘Lauranne’ (Goonetilleke et al., 2023). Crack-out percentage and globosity are inversely correlated, and a QTL for both traits has been identified in the same region of LG2. This correlation between these two traits has not been previously reported, and further results are needed to confirm whether this is a common phenomenon in almond or specific to this population. For the chemical kernel composition, it is noteworthy that there is no correlation between total fat content and the composition of fatty acids, but there is a high correlation between most of the different fatty acids, as previously observed (Font i Forcada et al., 2012). Additionally, there is a notable negative correlation between fat and protein content. The highest correlation between physical and chemical traits was observed for the fat content, which was inversely correlated with globosity (-0.45), and for protein content, which is inversely correlated with crack-out percentage (-0.39). These correlations are particularly significant, given that the only detected QTL for these traits was mapped to the same region of chromosome 2. Other interesting correlations of chemical traits have been explored, such as the relationship between oleic acid content and shelf life recently pointed out in rapeseed (Spasibionek et al., 2020) and almond (Sideli et al., 2024).

Several QTLs associated with kernel traits were identified by our study, some of which co-localize QTLs previously identified by other authors. For instance, the crack-out percentage QTL (qP-Crack2.1), located on chromosome 2, has been reported by previous studies for crack-out and shell hardness (Goonetilleke et al., 2018; Pavan et al., 2021; Sánchez-Pérez et al., 2007; Sideli et al., 2023; Pérez de los Cobos et al., 2023). Both ‘Marcona’ and ‘Marinada’ exhibited hard shells, however wide variation was found in their F₁ progeny. This suggests the dominance of hard shells and the presence of alleles in heterozygosity in the parental cultivars. The additive effect for crack-out percentage was 6.2% in ‘Marinada’ and 7.3% in ‘Marcona’. Meanwhile, in the MC × MI map, the difference between the highest and lowest mean values was 15.7%, closely matching the combined effect of both parents and supporting this hypothesis. Similarly, for kernel weight, a major QTL was identified at the beginning of LG4, the same chromosome where a previous QTL for this trait was identified in another F₁ population derived from the cross ‘R1000’ × ‘Desmayo Largueta’ (Sánchez-Pérez et al., 2007) but at a slightly different genomic position (approx. 3 Mbp apart). QTLs for kernel weight have been identified in various genomic regions across studies (Fernández i Martí et al., 2013; Goonetilleke et al., 2023; Pérez de los Cobos et al., 2023), suggesting that kernel weight is a quantitative trait controlled by several genes in the almond germplasm. This was corroborated in our study, where in addition to the main kernel weight QTL (qP-KWe4.1), an additional QTL, detected only in the ‘Marcona’ map, indicated the presence of multiple alleles associated with weight variation in these cultivars. Additionally, a QTL for kernel width (qP-KWidth_SA4.1) was also identified, overlapping with qP-We4.1, suggesting that width may influence the observed differences in fruit weight. These interactions and the numerous QTLs reported for this trait highlight the challenges in implementing effective breeding strategies for weight gain in breeding programs.

QTLs for various kernel morphology traits, including length, thickness, width, roundness, symmetry, and shoulder angle were identified at various genomic locations. Notably, major QTLs for kernel length and width were found 25 Mbps apart on chromosome 1, suggesting that different genetic mechanisms control these traits. qP-KLen1.1, associated with kernel length, had previously been mapped in a panel of 98 almond cultivars (Font i Forcada et al., 2015a) and in an F₁ population (‘Vivot’ × ‘Blanquerna’) (Fernández i Martí et al., 2013). Additionally, Goonetilleke et al. (2023) reported two overlapping QTLs for kernel width and length on LG1 of the ‘Lauranne’ map. However, this region is situated between the kernel length and width QTLs identified in the MC × MI population, suggesting that qP-KWidth1.1 is a new QTL associated with kernel width found in almonds. For roundness and symmetry, QTLs for the same traits were recently reported in the same genetic region of chromosome 1 (Goonetilleke et al., 2023). A QTL for kernel thickness was also identified on LG4 (qP-KThick4.1). Finally, another QTL for roundness was mapped on LG7 in a different position where a QTL for sphericity was reported on the same LG in a ‘Nonpareil’ × ‘Lauranne’ population (Goonetilleke et al., 2023).

In this population, we studied the segregation of the shoulder character, which is being analyzed for the first time in this study. We observed a strong correlation between kernel shoulder, width, and thickness, with a negative correlation with kernel length. Similar to the pattern seen with shape QTLs, QTLs associated with kernel shoulder exhibit a broader confidence interval compared to those for length or width, and similar to shape–related traits. This suggests that the QTL for kernel shoulder is likely influenced by both kernel length and width QTLs. To further validate QTLs for kernel length and width, it would be preferable to use a population where the kernel shoulder is not segregating.

Manual measurements for kernel shape traits and image analysis using artificial intelligence were employed in this study for QTL mapping. Despite the high correlation between both sets of phenotypic data, some QTLs identified through image analysis were not detected with manual measurements. This highlights the efficiency of using this new phenotyping method, not only in reducing data collection efforts but also in increasing the identification of new QTLs.

A cluster of overlapping QTLs located at the lower end of LG1 was identified for the fatty acids. These results indicate that a single molecular marker could predict the content of all these fatty acids. This is particularly significant given the prevalence of these acids in almonds, their role in health benefits associated with almond consumption, and their impact on kernel rancidity, which compromises flavor and postharvest shelf life (Barreca et al., 2020; Becerra-Tomás et al., 2019; Flankin and Mitchell, 2019). Although, these fatty acids were quantified in only one year; the fact that all acids mapped to the same region of LG1, along with the high percentage of variance explained by these QTLs (over 25% for most of the acids), underscores the reliability of the results despite validation with additional year data. Moreover, this region coincides with a QTL previously documented in almonds associated with fatty acids content and located at the bottom of LG1 in a biparental population (‘Vivot’ × ‘Blanquerna’) and a panel of 98 almond accessions (Font i Forcada et al., 2012, 2015b; Sideli et al., 2024). Additionally, this region of LG1 encompasses markers recently associated with different degrees of rancidity (Sideli et al., 2024). It is known that fatty acids and phenolic compounds are related to kernel rancidity, so the overlap of these two regions highlights the impact of these acids on controlling rancidity levels. Furthermore, the additive effect observed with these acids (up to 28.7 mg/g kernel weight for oleic acid) underscores the significant potential for improvement through breeding, which could result in new cultivars with an enhanced fatty acid profile.

We identified four QTLs associated with protein content on LGs 2, 3, 5, and 7. Among these, only qP-Protein3.1 and qP-Protein7.1 co-localize with QTLs mapped by previous studies (Font i Forcada et al., 2012, 2015b). However, the primary QTL for protein content discovered in this study, qP-Protein2.1, did not align with any previously identified QTL region, suggesting that the regulation of protein content involves multiple genes. For fiber and fat content, QTLs related to these traits are reported for the first time by our study. Both qP-Fiber3.1 and qP-Fat3.1 showed low significance, which could be related to their dependence on environmental conditions, as different stresses like drought and temperature have been linked to interannual variation in these compounds (Kodad et al., 2018). This is consistent with the low proportion of phenotypic variability explained by cultivars for fiber (15.7%) and fat (28.4%) contents in previous studies (Romero et al., 2011), highlighting the low heritability of these compounds.