Svs v8

Genetic diversity values were calculated by a neighbor-joining algorithm using TASSEL 5 (www.maizegenetics.net). The EIGENSTRAT algorithm [40 (link)] with the SNP and Variation Suite (SVS v8.8.5; Golden Helix, Bozeman, MT, USA, www.goldenhelix.com) was used for Principal Component Analysis (PCA). Observed nucleotide diversity (π) for various chromosomes was estimated with sliding-window analysis by using TASSEL v5.0 [41 (link)]. The fixation index (F_ST) was calculated by using Wright’s F statistic [42 ] with SVS v8.8.5 (Golden Helix, Bozeman, MT, USA, www.goldenhelix.com).

For GBS data, we considered only SNPs successfully mapped to the whole-genome sequence draft, because knowing the physical location of SNPs helps prevent spurious LD and, thereby, calling unreliable haplotype blocks. Mapped SNPs were further filtered by call rate >90%. Before studying LD decay, haplotype blocks were calculated for all markers by using the default settings in SVS v8.8.5. (Golden Helix, Inc., Bozeman, MT, USA, www.goldenhelix.com). Adjacent and pairwise measurements of LD for GBS data were calculated separately for SNPs in each scaffold. All LD plots and LD measurements and haplotype frequency calculations involved using SVS v8.8.5 (Golden Helix, Inc., Bozeman, MT, USA, www.goldenhelix.com).

Expected nucleotide diversity (π) for various chromosomes were estimated with sliding-window analysis by using TASSEL v5.0 as described [49 (link)]. Estimation of fixation index (F_ST) was based on Wright’s F statistic [50 (link)] with use of SVS v8.8.5 (Golden Helix, Inc., Bozeman, MT, USA, www.goldenhelix.com).

ROH detection analysis was carried out with a subset of SNPs, selected for having MAF > 0.15. The algorithm implemented in SVS v.8.8.3 (Golden Helix Inc.) was used to identify completely homozygous genomic stretches on chromosomes 1–8 with at least 15 SNP loci and with a minimal length of 100 Kb. As a measure of inbreeding, the ROH count and the ROH total length were computed for each individual. The ggplot2 R package⁵⁸ was used to visualize the distribution of the percentage of missing data per ROH and to perform a regression analysis between mean read depth per cultivar and ROH count per cultivar.
Individual inbreeding was also estimated using the F_PLINK inbreeding coefficient, which was computed using the LD-pruned marker dataset as input. Regression analysis between ROH count, or ROH total length, and F_PLINK coefficient, was performed using the ggplot2 R package⁵⁸ .

Genetic structure was studied using SNPs in approximate linkage equilibrium, which were obtained using the LD pruning algorithm in PLINK v.1.90p^{17 (link)}. This calculates pairwise R² for all marker pairs in sliding windows with a size of 50 markers and an increment of 5 markers and removes the first marker of pairs, in which R² < 0.5.
Analysis with the ADMIXTURE parametric model^{35 (link)} was performed with a number of ancestral populations (K) ranging from 1 to 15. One thousand bootstrap replicates were run to estimate parameter standard errors. The most suitable number of K was selected in correspondence with the lowest cross‐validation (CV) error. Cultivars were assigned to one specific ancestral population when the membership coefficient q_i for that cluster was >0.6. If not, they were considered admixed.
PCA on SNP data was performed using SVS v.8.8.3 (Golden Helix Inc.), and a three-dimensional plot was obtained using the top three components identified with default parameters of the additive model.
TreeMix (v1.12)^{36 (link)} was used to infer splits and mixtures among Italian, French, Spanish, and U.S. germplasm, testing a model with no migration, and models with all the three possible migration events among the four populations. The “get_f()” R function was used to obtain the variance explained by each model.

HM was carried out using options available in SVS v.8.8.3 (Golden Helix Inc.). Clusters of ROHs, defined as genomic regions of at least 100 Kb whose loci occur in ROHs of at least ten cultivars characterized at the phenotypic level, were identified. Repeated binary spectral clustering^{62 (link)} was used to trim boundaries of clusters of ROHs, in order to define homozygous regions highly overlapping among cultivars. Finally, a linear regression model was fit between clusters of ROHs and BLUPs, using the top five principal components as covariates to correct for population structure. The FDR correction was used to account for multiple testing and suggest an association for P < 0.1. BLUP means of cultivars either contributing or not contributing to clusters of ROHs associated with phenotypic traits were computed and compared using a two-tail Student’s t test. Genes included in clusters of ROHs identified by HM were retrieved by the Prunus dulcis (cv. Lauranne) v1.0 genome annotation available at the genomic database of Rosaceae^{28 (link),63 (link)}.