A total of 92 accessions of two Spinacia spp. was used in this study. Three accessions of S. turkestanica originating from Turkmenistan, another S. turkestanica accession of unknown origin, and 16 S. oleracea accessions were received from the National Plant Germplasm System (NPGS) of the United States Department of Agriculture (USDA). In addition, 72 S. turkestanica accessions collected in Uzbekistan and Tajikistan were obtained from the CGN, Wageningen University and Research (WUR). CGN collected the material under the Standard Material Transfer Agreement (SMTA) of the ITPGRFA. Information on each accession, including the country of origin, collection site (latitude, longitude, and altitude), and collection year is provided in Supplementary Table S1. For each accession, one randomly selected plant (representing an accession) was used for the molecular analyses.

Spinach plants were grown in RediEarth propagation mix (Sunshine Horticulture, Agawam, MA) in a greenhouse at the Washington State University (WSU) Mount Vernon Northwestern Washington Research and Extension Center, with the air temperature set at 22–24°C during the day and 18–20°C by night, and supplemental lighting provided for 10 h/day. Plants were fertigated daily with General Purpose Fertilizer 20-20-20 (Plant Marvel, Chicago, IL) injected into the irrigation water at a 1:100 ratio for applying a final nitrogen concentration of 200 ppm at each irrigation. A single leaf was harvested from each plant 35 days after planting, frozen immediately in liquid nitrogen, and stored at −80°C.

Genomic DNA was isolated from the single frozen leaf of each plant using the cetyl trimethylammonium bromide (CTAB) method (Murray and Thompson, 1980; Porebski et al., 1997), after the leaf was ground in liquid nitrogen with a mortar and pestle. The DNA quality was checked on a 1% agarose gel, quantified using a Qubit, and submitted to the University of Wisconsin Madison Biotechnology Center (https://www.biotech.wisc.edu/) for sequencing, where DNA quality and integrity were re-evaluated using the Quant-IT PicoGreen fluorescent dye (Thermo Fisher, Waltham, MA). The GBS method of Elshire et al. (2011) was used to sequence the samples after digesting genomic DNA with the ApeKI restriction enzyme, as described by Bhattarai et al. (2020). Digested DNA fragments were ligated with unique barcodes and Illumina adapters, and the samples were pooled in equal proportion to construct GBS libraries, as described by Elshire et al. (2011). Finally, the 96-plex GBS libraries were amplified, purified, and sequenced as 150 bp paired-end reads on an Illumina NovaSeq machine (Illumina, San Diego, CA).

Using the Skewer program (Jiang et al., 2014), the reads were pre-processed to remove sequencing adapters and to filter out low-quality bases for a minimum quality of Q20. Filtered, good quality reads were de-multiplexed and aligned to the spinach reference genome (Xu et al., 2017; http://www.spianchbase.org) using Bowtie two software (Langmead and Salzberg, 2012). The TASSEL GBS v2 pipeline was used to remove barcodes, filter for quality, and call SNPs (Bradbury et al., 2007; Glaubitz et al., 2014). The SNPs were filtered further using VCFtools v0.1.15 (Danecek et al., 2011) to remove indels, remove minor alleles with frequency (MAF) < 0.05, retain only bi-allelic SNPs, achieve a minimum SNP quality score (minQ) < 20 and a minimum genotype read depth (minDP) of 10, remove SNPs missing from >20% of the accessions, and remove plants (accessions) with >20% missing data. Final filtration was used to exclude SNPs with >20% missing data, a minDP of 15, and a MAF <5%. The filtered SNPS were pruned for linkage disequilibrium (LD) in PLINK v1.9 (Chang et al., 2015) using the indep-pairwise 50 5 0.2 option to remove correlated pairs of SNPs. The SNPs were separated into two files, one for each of the S. oleracea and S. turkestanica species. Both files were checked for common SNPs, and the common SNP was genotyped across all accessions for both species. Final filtered SNP distribution across the six spinach chromosomes was determined using the CMplot package in R. The filtered datasets were then used for diversity, phylogenetic, structure, and selective sweep analyses.

The structure of accessions of the two Spinacia spp. was analyzed using the software STRUCTURE 2.3.4, with individual accessions assigned to genetic clusters, hereafter called groups, based on inferred genetic ancestry (Pritchard et al., 2000). Structure implements a Bayesian model-based clustering method which assigns multi-locus accessions to a number of user-defined groups (K), and is based on maximization of linkage equilibrium (LE) within groups but minimization of LE among groups. The structure analysis parameters were set to an admixture model, with K ranging from 1 to 10 using five iterations, a burn in period of 100,000, and a Markov Chain Monte Carlo (MCMC) run length of 100,000. The K was determined by considering various factors suggested by Camus-Kulandaivelu et al. (2006). First, Structure Harvester version 0.6.4 was used to determine K, as described by Evanno et al. (2005) [http://taylor0.biology.ucla.edu/structureHarvester/]. The resulting proportion of membership coefficients (Q matrices) for each accession was used to draw a bar plot to visualize clustering of the spinach accessions. True groups were identified as the maximum value of ΔK, based on the rate of change of the natural log probability of the data. Spinacia spp. were assigned to individual groups (Q) based on an assignment value of 75%. Second, the inferred ancestry assessments obtained from Structure were judged based on prior knowledge of the geographic location of the collection site of each accession of both Spinacia spp., provided by the respective gene banks. A total of 7,065 SNPs common to the 76 S. turkestanica and 16 S. oleracea accessions was used to analyze the population structure of the accessions of these two species. To understand better the group structure and clustering pattern of these accessions, principal component analysis (PCA) was completed using PLINK v1.9 (Chang et al., 2015), and plotted in R. Genetic relationships among the accessions were inferred based on a neighbor-joining phylogenetic analysis in MEGA7 (Kumar et al., 2016) with 200 bootstraps.

The SNPs were used to determine the genetic diversity and genetic structure of the accessions of the two Spinacia species. The summary statistics of genetic diversity (GD) and polymorphic information content (PIC) of SNPs were calculated using PowerMarker software V 3.25. (Liu and Muse, 2005). The PIC of SNPs was calculated using the following formula, according to Botstein et al. (1980): P I C = 1 −   ∑ j = 1 n P i j 2 − ∑ j = 1 n = 1 ∑ k = j + 1 n 2 P i j 2 P i k 2 where, P _ij and P _ik are the frequencies of the jth and kth alleles, respectively, of bi-allelic SNP marker i.

The number of different alleles (Na), number of effective alleles (Ne), Nei’s genetic diversity (h) (Nei, 1973; Nei, 1978), unbiased genetic diversity (uh), and Shannon’s information index (I) were calculated using GenAlEx 6.3 (Peakall and Smouse, 2006). Analysis of molecular variance (AMOVA), group genetic differentiation (Fst), Nei’s unbiased genetic distance (D), and Nei’s unbiased genetic identity (I _d ) were also estimated using GenAlEx 6.3. Gene flow (Nm) was calculated using the formula Nm = [(1/Fst)-1]/4, according to Peakall and Smouse (2006). The AMOVA (Excoffier et al., 1992) was used for hierarchical partitioning of genetic variation among the groups, and among individuals within the groups of the two Spinacia spp., using G-statistics. In total, 999 permutations determined the fixation indices (Fst) at a level of significance of p < 0.001. Nucleotide diversity (π) was calculated using VCFtools v0.1.15 (Danecek et al., 2011) with a 1-kb window size across the genome (window-pi 1,000).

Linkage disequilibrium (LD) was calculated using SNP pairs within a 200 Kb window, and plotted using PopLDdecay3.4.1 (Zhang et al., 2019) for each S. oleracea and S. turkestanica, and for the Q2 and Q3 groups of S. turkestanica (see details in the Results section). Again, for LD analysis between the Q2 and Q3 groups of S. turkestanica, accessions with <0.75 membership assignment (admixed between Q2 and Q3) of the group structure were excluded. As a result, of 76 S. turkestanica accessions, only 63 (33 from Q2 as the query panel, and 30 from Q3 as the reference panel) were used. The physical distance of the LD curve intersecting at the critical value of r ² = 0.20 was used as LD decay for each species and for groups within S. turkestanica. Screening of the spinach genome for selective sweep regions was first performed by comparing allele frequency differentiation between the S. turkestanica and S. oleracea accessions, following the method described by Xu et al. (2017). Then, since the population structure analysis described above revealed two distinct groups, Q2 and Q3, among the S. turkestanica accessions, screening for selective sweep regions associated with the Q2 and Q3 groups of S. turkestanica was carried out by modelling the likelihood of multi-locus allele frequency differentiation between the two groups using XP-CLR v1.0 (Chen et al., 2010). XP-CLR detects selective sweep regions by modeling the likelihood of multi-locus allele frequency differentiation between two groups. The XP-CLR was run using xpclr v1.1.2 (https://github.com/hardingnj/xpclr) for each pseudochromosome with 50 Kb sliding window, 10 Kb step size, and setting the maximum number of SNPs in each window to 50 (--size 50,000 --step 10,000 --maxsnps 50). The adjacent windows (<10 Kb) with high XP-CLR scores (top 1%) were grouped into a single region representing a single selective sweep region. The candidate genes within the selective sweep regions were also identified. XP-CLR analysis was performed using the 76 S. turkestanica accessions as a reference panel and 16 S. oleracea accessions as a query panel to identify selective sweeps. Of the 63 S. turkestanica accessions that remained after filtering for admixture, 33 accessions in Q2 and 30 in Q3 were used as the reference panel in the XP-CLR analysis. Similarly, given the evidence (see Results) that the S. oleracea accessions (Q1) are most closely related to the Q2 group of S. turkestanica, an analysis of selective sweeps was also computed between the 16 Q1 accessions as the query group and 33 Q2 accessions as the reference group.

Approximately 331.2 million raw reads were generated from the Illumina NovaSeq run for the 92 Spinacia accessions. After filtering for sequencing adapters, low-quality bases, and de-multiplexing to extract reads matching the sample barcodes, approximately 300.7 million reads were retained with an average of 3.1 million and a median of 3.1 million reads per accession. Using the TASSEL GBS v2 pipeline, 210,968 SNPs were identified across the six chromosomes. After the SNPs were filtered for indels and a MAF <0.05, 105,328 SNPs were retained. Further filtering to keep only bi-allelic SNPs resulted in retention of 103,704 SNPs.

When the SNP dataset was separated into two files, one for each species, with filtering for 20% missing SNP data, 9,456 and 7,629 SNPs were retained for S. oleracea and S. turkestanica, respectively. Across both species, there were 20.7% missing calls, and the rest of the SNPs were retained for genetic analysis. Among the filtered SNPs, 7,065 SNPs common to all 92 Spinacia accessions that met the filtration criteria were retained. The distribution and density of these SNPs across the six chromosomes of S. turkestanica and S. oleracea are presented in Figure 1, with 832, 839, 1,915, 1,748, 906, and 825 SNPs located on chromosomes 1 through 6, respectively. The SNP density ranged from 0 to >40 SNPs/Mb physical distance. These SNPs were well distributed across all six chromosomes (Figure 1).

The genetic diversity of the 7,065 SNPs common to the 16 S. oleracea and 76 S. turkestanica accessions ranged from 0.10 (40 SNPs) to 0.50 (3,418 SNPs), with an average of 0.35 (Supplementary Figure S1A). The genetic diversity of 79% of these SNPs was >0.3. The PIC of these SNPs ranged from 0.10 (160 SNPs) to 0.40 (3,798 SNPs), with an average of 0.28 for the S. oleracea and S. turkestanica accessions (Supplementary Figure S1B). The PIC of 75% of these SNPs was >0.3.

The AMOVA of the 76 S. turkestanica and 16 S. oleracea accessions with 7,065 common SNPs showed 0.5% estimated variance accounted for by the two species, while 99.5% was accounted for by individual accessions within each species (Table 1). The average Fst and gene flow (Nm) was 44.5 among three groups of the two species. The greatest Fst, 0.012, was between groups Q2 and Q3, both of which consisted of accessions of S. turkestanica, followed by 0.011 for group Q1 of S. oleracea and Q3 of S. turkestanica, and then 0.008 between the Q1 and Q2 groups (Table 2). The genetic distance, D, followed the same pattern, with the largest D between groups Q2 and Q3, and the smallest between groups Q1 and Q2. The genetic identity, I _d , was greatest for Q1 and Q2 (0.993), followed by Q1 and Q3 (0.989), and Q2 and Q3 (0.988). Consequently, the highest Nm was between Q1 and Q2 (30.13), followed by Q1 and Q3 (21.83), and Q2 and Q3 (19.77) (Table 2). The mean number of different alleles (Na), number of effective alleles (Ne), diversity index (h), unbiased diversity index (uh), Shannon’s Information Index (I) and nucleotide diversity (π) ranged from 1.988 to 1.997, 1.607 to 1.616, 0.347 to 0.351, 0.355 to 0.357, 0.516 to 0.523, and 0.9223 × 10⁻³ to 0.9371 × 10⁻³ respectively (Table 3). The mean nucleotide diversity (π) of the 92 Spinacia accessions was estimated to be 0.8998 × 10⁻³. Overall, of the three groups, Q2 was the most diverse (h = 0.351, uh = 0.356, and I = 0.523), followed by Q3 (h = 0.349, uh = 0.355, and I = 0.520) and Q1 (h = 0.347, uh = 0.357, and I = 0.516). However, the unbiased diversity index (uh) of Q1, which largely comprised S. oleracea accessions, was greater than that of the Q2 and Q3 groups (Table 3).

The STRUCTURE analysis based on 7,065 SNPs common to the S. oleracea and S. turkestanica accessions revealed three main groups, Q1, Q2, and Q3 (Figure 2), that comprised 18, 43, and 31 accessions, respectively (Figure 3A). The Q1 group largely comprised cultivated spinach, S. oleracea, and admixtures of both species, while the Q2 and Q3 groups comprised S. turkestanica accessions only. Multivariate analyses using UPGMA and PCA supported results of the structure analyses, clustering a large majority of the accessions of Spinacia into three groups, Q1 (11 accessions), Q2 (34 accessions), and Q3 (30 accessions), with another 17 accessions admixed (Figures 3B,C). The accessions AM45 (United States), AM239 (the Netherlands), AM268 (Macedonia), AM277 (United Kingdom), AM316 (Denmark), AM330 (Nepal), AM360 (United States), and Viroflay (France) had >0.970 inferred membership in the S. oleracea group, while another eight S. oleracea accessions had varying degrees of admixture with S. turkestanica (Figure 3A).

The distribution of spinach accessions based on their original geographic location of collection is presented in Figure 4, with detailed information on the two Spinacia spp. provided in Supplementary Table S1. The Q1 group was dominated by S. oleracea landraces (accessions) collected from Afghanistan, China, Iran, India, Nepal, and Pakistan, as well as nine S. turkestanica accessions admixed with S. oleracea accessions (Figure 4). Of these nine admixed S. turkestanica accessions in Q1, five (St22, St36, US-St4, US-St5, and US-St8) were collected from Turkmenistan. Interestingly, with few exceptions, group Q2 comprised S. turkestanica accessions originating from Uzbekistan, while a majority of group Q3 comprised S. turkestanica accessions originating from Tajikistan (Figure 4). The exceptions included three accessions (St39, St29, and St58) originating from the southeastern part of Uzbekistan, which had varying degrees of admixture of groups Q2 and Q3.

The LD plot showed decreasing LD between markers with the increase in physical distance on the chromosomes (Figure 5). The LD decay plot revealed a rapid rate of LD decay in both Spinacia species. The LD decay was around 9 Kb in S. turkestanica and around 12 Kb in S. oleracea at r ² = 0.2 (Figure 5A). The S. turkestanica groups Q2 and Q3 did not show such differences in LD decay rate (Figure 5B). The analysis of selective sweeps of the S. oleracea and S. turkestanica genomes, determined by XP-CLR analysis, revealed a total of 20 regions, including 3, 2, 9, 1, 3, and two regions in chromosomes 1, 2, 3, 4, 5, and 6, respectively (Table 4; Supplementary Table S2). The highest XP-CLR score was for the 30–32 Mb region of chromosome 1, with 35 SNPs under selection (12 SNPs in chromosome region 1.3 with an XP-CLR score of 24.75009, and 23 SNPs in region 1.2 with an XP-CLR score of 0.00597). Another important selective sweep region was at 98 Mb on chromosome 3, where 23 SNPs are under selection (12 SNPs in region 3.8 with an XP-CLR score of 18.20349, and 11 SNPs in region 3.9 with a score of 11.40893).

Several of the selective sweep regions of Q1 vs. Q2 accessions were common to the Q1 vs. all S. turkestanica accessions (Q2 and Q3), including sweep regions 1.3, 2.1, 3.2, 3.8, 3.9, 4.1, 5.1, 5.2, 6.1 and 6.2 (Table 4). The analysis of selective sweeps of the two groups of S. turkestanica, Q2 vs. Q3, showed a total of seven regions spread over chromosomes 2, 3, 5, and 6 that played important roles in the differentiation of these two groups. Of these selective sweeps, 10 SNPs at the 52 Mb region of chromosome three were common to both S. oleracea and S. turkestanica and had an XP-CLR score of 19.99774. The next highest XP-CLR score, 5.23574, was for the selective sweep region at 10 Mb on chromosome 5, which encompassed 12 SNPs (Table 4).