Plant materials consisted of 155 accessions of L. sativa and 5 of the closely related wild species L. serriola, which were part of the germplasm panel established in the frame of the ECPGR European Evaluation (EVA) Lettuce Network (ECPGR, 2023). Plant materials originated from the germplasm collections of four institutions: the Institute for Plant Genetic Resources “K.Malkov” (Sadovo, Plovdiv district, Bulgaria), the Centre for Genetic Resources, the Netherlands (CGN, Wageningen, Netherlands), the Unité de Génétique et Amélioration des Fruits et Légumes, Plant Biology and Breeding, INRAe (GAFL, Avignon, Montfavet Cedex, France), and the Nordic Genetic Resource Center (Nordgen, Alnarp, Sweden). Genotypes encompassed cultivars, breeding materials, and landraces originating from a total of 10 different countries in Europe, America, and Asia. Different horticultural types were represented ( Figure 1 ), including Butterhead (54), Iceberg (46), Cos or Romaine (17), Batavia or Summer/French Crisp (11), Crisp (10), Loose leaf (9), Oak leaf (4), Latin (3), and Lollo (1), as well as wild L. serriola (5) (also known as prickly lettuce). A detailed list with all available information on the assayed accessions is provided in Supplementary Table 1 .

For probe design, a dataset including whole-genome resequencing data of 131 L. sativa accessions (Wei et al., 2021) was considered ( Supplementary Table 2 ). Raw sequence data were retrieved from the FTP site of the China National Gene Bank Sequence Archive (CNSA) repository (Guo et al., 2020). Variants (SNP and INDEL separately) were selected by filtering those present in the dataset with a minimum allele count of 3 (i.e., one homozygous accession and one heterozygous or three heterozygous accessions). The lettuce reference genome (L. sativa cv Salinas V8) and its annotation were retrieved from https://lgr.genomecenter.ucdavis.edu/Home.php and all gene coordinates were extended by 5,000 bp upstream and 1,000 bp downstream. All selected genomic variants were intersected with these gene coordinates and labeled as gene-space variants. A panel of 50k target sites was then built by imposing a minimum distance of 3,000 bp for variants on the gene-space and a minimum distance of 200,000 bp in the intergenic regions. After two rounds of design, a final panel of 41,547 targets were successfully identified by unique probes. Each probe consisted of a 40-bp sequence. SNP calling was enabled 460 bp downstream of the probe.

Genomic DNA was isolated from young leaves of a single individual per accession using a NucleoSpin Plant II Mini kit (Macherey-Nagel GmbH & Co. KG., Düren, Germany. DNA concentration was measured using the Qubit 2.0 Fluorometer (Thermo Fisher Scientific, Waltham, MA, USA). Libraries were prepared using the “Allegro Targeted Genotyping” protocol from NuGEN Technologies (San Carlos, CA), using 10 ng/μl of DNA as input and following the manufacturer’s instructions. Libraries were quantified using the Qubit 2.0 Fluorometer, and their size was checked using the High-Sensitivity DNA assay from Bioanalyzer (Agilent technologies, Santa Clara, CA) or the High-Sensitivity DNA assay from Caliper LabChip GX (Caliper Life Sciences, Alameda CA). Libraries were quantified through qPCR using the CFX96 Touch Real-Time PCR Detection System (Bio-Rad Laboratories, Hercules, CA) and sequenced on the Illumina NovaSeq 6000 (Illumina, San Carlos, CA).

Demultiplexing of raw sequencing data and base calling (BCL files into FASTQ files) were performed with the Illumina bcl2fastq2 Conversion Software v2.20 (Illumina, San Carlos, CA). Read quality check and adapter trimming were carried out using ERNE v1.4.6 (Del Fabbro et al., 2013) and Cutadapt (Martin, 2011), both with default parameters. Alignment to the reference genome L. sativa cv Salinas V8 (Reyes-Chin-Wo et al., 2017) was done using the Burrows–Wheeler Aligner BWA-MEM v0.7.17 (Li and Durbin, 2009) with default parameters and selection of uniquely aligned reads (i.e., reads with a mapping quality >10). SNP calling was obtained using gatk-4.0 (DePristo et al., 2011) following the software best practices for germline short variant discovery. SNP calling was limited to the regions (460 bp) that were previously defined as downstream of each enrichment probe.

All analyses were implemented in GATK Best Practices v4.1.2.0 (Van der Auwera and O'Connor, 2020) and included the following steps: (i) per-sample variants calling on target regions using HaplotypeCaller with default parameters to create a GVCFs file for each sample; (ii) GVCFs consolidation across multiple samples using GenomicsDBImport with default parameters and target intervals in order to improve scalability and speed for further joint genotyping; (iii) joint genotyping using GenotypeGVCFs with default parameters to produce a set of joint-called variants; (iv) Selection of SNPs using SelectVariants and quality filtering of SNPs using VariantFiltration (filter expression used: QD< 2.0 || MQ< 40.0 || MQRankSum< −12.5). A 1,911,467 biallelic SNPs matrix was obtained. The extra filtration of the VCF was performed with bcftools by setting all data points with fewer than five reads in coverage to a missing data genotype (./.) and retaining only records where a minimum of 96 samples reported a coverage above 10 reads. In total, 835,426 SNPs were obtained. For downstream analysis, 81,531 SNP sites were retained with minor allele count = 3, max missing 0.5, minQ = 30, and minor allele frequency = 5%. VCFtools version 0.1.17 (Danecek et al., 2011) was used. Pattern of nucleotide diversity (p) was estimated in non-overlapping sliding windows with a size of 1 kbp in VCFtools. Functional annotation of the identified variants associated genes was performed using SnpEff (version 3.1) (Cingolani, 2022).

Genetic diversity summary of the SNP matrix was performed by the Geno Summary tool implemented in Tassel v5.2.15 (Bradbury et al., 2007). Considering the biallelic nature of SNPs, expected heterozygosity according to Hardy–Weinberg equilibrium (H) was calculated according to the formula

where p and q each represent the frequency of the different alleles for each SNP.

The polymorphic information content (PIC) was calculated according to the formula (Shete et al., 2000)

Population structure was determined using the model-based ancestry estimation obtained with ADMIXTURE software (Alexander et al., 2015) with K ranging from 1 to 15. One thousand bootstrap replicates were run to estimate parameter standard errors. Tenfold cross-validation (CV) procedure with five iterations was performed, and CV scores were used to determine the best K value. Individuals were considered to belong to a specific K population if its membership coefficient (qi) was ≥0.5, whereas the genotypes with qi lower than 0.5 at each assigned K were considered as admixed. A neighbor-joining phylogenetic tree was built using the Jones–Taylor–Thornton (JTT) model with 1,000 bootstraps. Analyses were conducted in MEGA X software (Kumar et al., 2018). Principal component analysis (PCA) was performed in Tassel v5.2.15 and the biplot was drawn using the ggplot2 R package (Wickham, 2016).

The phenotypic traits were surveyed across five locations (Eyragues, Avignon, France; La Ménitré, France; Les Evouettes, Port-Valais, Switzerland; Rheinau, Switzerland; and Thessaloniki, Greece) during the 2020–2022 spring seasons. Plants were grown in a randomized block design with three replicates. Field trials were conducted using the standard agricultural practices for the local area of cultivation. Four traits were assayed including (i) seed color (1 = white/cream, 2 = yellow, 3 = brown, 4 = black), (ii) outer leaf color before bolting stage (1 = yellow green, 2 = green, 3 = gray green, 4 = blue green, 5 = red green), (iii) leaf anthocyanin content before bolting stage (0 = absent, 3 = weak, 5 = medium, 7 = strong), and (iv) bolting time (number of days from sowing to bolting).

Genome-wide association analysis was performed in 155 L. sativa genotypes. Six models were used including the general linear model (GLM) (Loley et al., 2013), the mixed linear model (MLM) (Zhang et al., 2010), the multi-locus mixed linear model (MLMM) (Segura et al., 2012), the compressed mixed linear model (CMLM) with population parameters previously defined (P3D) (Zhang et al. in 2010), the fixed and random model circulating probability unification model (FarmCPU) (Liu et al., 2016), and the Bayesian-information and Linkage-disequilibrium Iteratively Nested Keyway model (BLINK) (Huang et al., 2019). All models included the population structure as a covariate. The kinship was estimated using the identity by state (IBS) for accounting relationships among individuals. Phenotypic data from independent experiments were implemented. The significance threshold for marker–trait association was determined after Bonferroni multiple test correction with genome-wide α = 0.05. Considering 81,531 SNPs, the marker was considered significant when the p-value was less than 6.212 (−log₁₀P = 6.133 × 10⁻⁷). GLM and CMLM were computed in Tassel v 5.2.82 (Bradbury et al., 2007). MLM, MLMM, FarmCPU, and BLINK were calculated with the GAPIT R package (Wang and Zhang, 2021). Manhattan and quantile–quantile (Q–Q) plots for GWAS results were produced using the R package CMplot. The chromosomal location of the genome-wide significantly associated SNPs was displayed using PhenoGram (https://ritchielab.org/software/phenogram). Significant association signals were checked for their physical position on the L. sativa (cv. Salinas) V8 genome. The information about predicted genes was downloaded from the Lettuce genome browser v8.0 (https://phytozome-next.jgi.doe.gov/jbrowse/). Underlying genes and their functions were determined according to Reyes-Chin-Wo et al. (2017).

Based on SNP data retrieved from 131 lettuce raw sequences (Wei et al., 2021) and on the alignment to reference genome L. sativa cv Salinas V8, 41,547 probes were designed, of which 1,707 (4.1%) were localized in intergenic regions and 39,840 (95.9%) within genes ( Supplementary Table 3 ). The average coverage of the total set of probes was 77.1×; for those located within intergenic regions, 93.4×; and for those within genes, 76.6× ( Supplementary Figure 1 ). The SPET panel showed an average distribution of one probe per 55.5 kilobase pair (kbp). Regarding inter-probe distance, 27% of the probes were more than 50 kbp apart, while the largest gap was 3.2 mega base pair (Mbp) on chromosome 3 ( Figure 2 ). The sequencing of SPET libraries in the 160 study samples produced a total of 668,695,867 paired end raw reads corresponding to an average of 4,179,349 read pairs per sample ranging from 2,000 to 21 million and a mean depth of 79.7× ( Supplementary Table 4 ). The mapping rate on the whole genome was on average 88%, and only eight samples had an average below 70% ( Supplementary Figure 2 ).

By applying stringent filtering criteria, we identified 81,531 SNPs ranging from 5,291 on chromosome 6 to 13,920 on chromosome 4 ( Table 1 ). SNPs were predominantly located within transcript regions, covering over 65% of the gene space in all chromosomes. SNP effect analysis showed that the majority of SNPs (88.08%) have a possible modifier effect, while the rest exhibited low (6.96%), moderate (4.82%), and high (0.14%) impacts ( Supplementary Table 5 ). Within gene space, SNPs were mostly localized in upstream and downstream gene regions (27.45% and 16.27%, respectively). SNPs in exons and introns were 11.37% and 7.47%, respectively ( Supplementary Figure 3 ). The average density corresponded to one SNP every 28.99 kbp across the nine chromosomes, ranging from 21.48 kbp on chromosome 1 to 36.35 kbp on Chr 6. Across the whole set, PIC values ranged from 0.033 to 0.375 (data not shown) with a mean of 0.240. The minimum average PIC value was encountered on Chr 6 (0.226), while the maximum value was found in Chr 2 (0.258). Chr 6 and Chr 2 exhibited the lowest and highest nucleotide diversity with 4.681^e-4 and 6.459^e-4, respectively. On average, heterozygosity was 0.292, reaching values above 0.300 only on chromosomes 2 and 5. The observed transitions/transversions ratio was 2.12 ( Supplementary Figure 4A ). In particular, among transition events, C > T and G > A were the most abundant (18.697% and 18.225%, respectively), whereas C > A and A > T abounded within transversion events (4.826% and 4.506%, respectively). The allele content of the SNP matrix was balanced, being on average represented for 70% by the four nucleotide bases in homozygosity state ( Supplementary Figure 4B ).

An admixture-based clustering model implemented in the software ADMIXTURE (Alexander et al., 2015) was used to infer the genetic structure of the studied germplasm. Using the entire SNP dataset, results of CV error suggested six different clusters ( Supplementary Figure 5 ) representing the most likely number of subpopulations (K) ( Figure 3 ).

The subpopulations reflected to some extent a differentiation based on cultivar typology rather than country of provenance ( Figure 3 ; Supplementary Table 6 ). The first cluster (K1) grouped 21 accessions, mostly Iceberg and Cos lettuce types from Bulgaria. Butterhead were mostly grouped in clusters 2 (K2) and 6 (K6) and represented 62% and 67% of the total individuals within each cluster, respectively. The subpopulation 3 (K3) included several Batavia and Crisp types whereas Oak leaf types were included in cluster 4 (K4) together with Iceberg and Loose leaf types. Among the different cultivar types, Iceberg accessions were clustered in several subpopulations. Lactuca serriola accessions were grouped separately from the rest in a distinct group (K = 5). Thirty-two accessions belonging to 8 out of the 10 considered cultivar types were classified as admixed, as they showed values for the highest cluster membership coefficient (qi) lower than 0.5. The Fixation Index (F_ST ) values, measuring the population (K) differentiation based on SNP data, are reported in Table 2 .

The highest F_ST values were found between K5 and the other subpopulations, thus confirming the differentiation of the wild L. serriola from the cultivated L. sativa. The lowest divergence was found between clusters 1 and 3 (F_ST = 0.265) mostly comprising the same type of cultivars. Considering the average q-value at K = 6 ( Figure 4 ), the analysis showed how among the most represented cultivars, iceberg types were included in five out of the six detected clusters while butterheads were included in clusters 2, 3, and 6. Batavia, Crisp, and Lollo as well as Loose and Oak leaf types were mostly represented by clusters 3 and 4, respectively. The average heterozygosity of the accessions was on average lower than 4% in all cultivated variety groups ( Figure 5 ). Prickly lettuce accessions showed an average heterozygosity of 4.69% with values ranging from 4.64% to 4.99%. The same trend was observed among different subpopulations based on admixture analysis (data not shown). Twelve accessions belonging to Butterhead (7) and Iceberg (5) exhibited heterozygosity higher than 5% with values up to 6.61% (Butterhead) and 10.08% (Iceberg). Only a single accession, representing the Cos horticultural type, showed a relatively high heterozygosity of 16.11%.

Phylogenetic clustering and PCA were performed to find patterns of genetic variation among accessions. The phylogenetic network using the neighbor-joining method was generally in agreement with Admixture analysis. Two main subpopulations were detected. Group I mostly included butterhead types ( Figure 6A ) from the clusters K2 and K6 ( Figure 6B ).

Group II consisted of several icebergs, loose leaf, and cos types from clusters K1, K3, and K4 ( Figures 6A, B ). All prickly lettuce (L. serriola) genotypes were grouped closely together according to the cluster K5. The distribution of the accessions in the PCA bi-plot graph corroborated population structure analysis highlighting, among L. sativa accessions, a clustering of butterhead types compared to the rest ( Figure 7A ). Prickly lettuce genotypes were grouped apart on the second component, thus confirming the observed subpopulation in the ancestry analysis. A slight differentiation between French and Bulgarian germplasm was observed. Interestingly, several close relationships were found between Italian and Bulgarian accessions. Although more admixtures were found when the geographical provenance was considered, a general differentiation was observed between germplasm retrieved from Western and Eastern Europe ( Figure 7B ).

Genome-wide association scans using six models detected a total of 306 significant SNP–trait associations (STA) ( Supplementary Table 7 ) distributed across all chromosomes except for chromosome 6. The majority of STA were detected for seed color and leaf anthocyanin content: 133 and 117, respectively. Fifty-eight percent of associations were identified with the GLM, whereas among the five multi-locus models used, MLM and CMLM highlighted the highest number of association signals. Only for bolting time was no association found with multivariate models, except FarmCPU. Considering all models, chromosomes 5, 7, and 9 held over 94% of the STA, showing further several colocalizations. Manhattan plots showing the associations, their chromosomal positions, and Bonferroni threshold are shown in Figure 8 , and Q–Q plots for multi-model GWAS and physical position of STA are shown in Figure 9 . Two main clusters were found for leaf anthocyanin content and outer leaf color in a 160-kb region at 86 Mbp on chromosome 5 as well as in a 2-Mbp region at 150–152 Mbp on chromosome 9. Furthermore, different significant SNPs were detected on chromosome 7 for seed color in a 3-Mbp region at 49–52 Mbp position and for bolting time at 164 Mbp.

In order to narrow down to potential GWAS hotspots, we considered the top-ranked SNPs within each model ( Table 3 ). For seed color, five out of the six models detected the strongest signal at 50.40 Mbp on chromosome 7 in an intergenic region at 19.55 kb to an Atp-dependent rna helicase DEAH5. The percentage of phenotypic variation explained (PVE%) by each locus ranged from 0.03% to 45.26%. Only with the CMLM model was the highest peak found 147 kb downstream to the previous one (chromosome 7, 50.54 Mbp) and in correspondence to Cytokinin dehydrogenase 3 . For leaf anthocyanin content, five models detected a robust association on chromosome 5 at 86.12 Mbp in correspondence to Phototropin-2 with a PVE% ranging from 4.55% to 15.26%. In addition, all models detected the strongest STA on chromosome 9 at 152.91 Mbp within an MLO like protein 11. Also, for outer leaf color, the strongest associations were in both chromosome 5 and 9, at ~27 kbp distance from those identified for leaf anthocyanin content. For leaf color, a Signal peptidase complex subunit 3B was the candidate gene identified with GLM, CMLM, and FarmCPU on chromosome 5 at 86.15 Mbp. The three models exhibited a PVE% ranging from 2.78 to 20.32. Furthermore, all models detected the strong STA at 152.88 Mbp on chromosome 9 in correspondence to a General transcription factor 3C polypeptide 6 , with a PVE% ranging from 10.46% to 47.65%.

For bolting time, only GLM, BLINK, and FarmCPU revealed the strongest STA in an intergenic region at 1.77 kbp from FAR1-related sequence 10 located on chromosome 7 at 164.43 Mbp. The three models exhibited a PVE% ranging from 18.71% to 48.65%.

In this work, we investigated the effectiveness of SPET as a tool for high-throughput genotyping in lettuce. This method has been developed recently, but so far, very little information about how well it performs in plants is reported. To that end, we developed and validated a novel SNP panel enriched of intraspecific SNPs from 131 resequenced genomes and consisting of over 40,000 probes across the lettuce genome. The potentialities of SPET rely on the high-efficiency enrichment of targeted loci and the high scalability of up to thousands of probes in a single reaction (Scaglione et al., 2019). In addition, it offers the possibility of discovering novel SNPs by sequencing the genomic regions surrounding the target SNPs. Compared to other genotyping strategies for reducing genome complexity, this method offers full control of target sites, thus broadening the investigation of variation within genomic regions with a functional role. Furthermore, the possibility to detect SNPs within probe-defined regions improves reproducibility, thus enabling one to implement and/or compare genomic information from different genotyping experiments. Our main goal was to determine the applicability of SPET for assessing the diversity of a heterogeneous germplasm collection of lettuce including genotypes belonging to different horticultural types with diverse geographic origins. This work was done as part of the ECPGR European Evaluation Network (EVA) with the goal of improving the knowledge of crop genetic diversity and exploiting it to breed more resilient crops that can meet the major problems facing agriculture in the upcoming years (FAO, 2021; ECPGR, 2023). A more efficient use of crop diversity is essential for genetic improvement, management, and conservation of germplasm resources. The sequenced dataset comprised an average of 4 million SNPs per sample, which has been indicated to be adequate for processing several thousands of probes (Scaglione et al., 2019). Compared to the 25K SPET panel reported in peach (Baccichet et al., 2022) and the 5K SPET panel described for tomato, eggplant, and oil palm (Barchi et al., 2019; Herrero et al., 2020), the 40K SPET assay designed in lettuce provides a higher number of SNPs covering up to 96% of gene-rich regions. Our findings demonstrated how effective SPET is compared to other genotyping techniques for detecting SNPs within coding regions. In fact, prior studies in lettuce utilizing genotyping by sequencing revealed that the proportion of SNP loci within genic regions ranged from 0.94% to 27.6% (Park et al., 2021; Park et al., 2022).

We detected over 80,000 high-quality polymorphisms that were analyzed to determine population ancestry, phylogenetic relationships, and principal components among the EVA lettuce accessions. The three approaches were complementary, thus supporting the interpretation of results. In agreement with earlier findings (Park et al., 2021; Park et al., 2022; Simko, 2009; Stoffel and van Leeuwen, 2012), no admixture was found between L. sativa and L. serriola. Indeed, the accessions of the two species were clearly separated. This evidence promotes the potentiality of SPET for phylogenetic studies, as already observed in aubergine and tomato (Barchi et al., 2019). Population structure and phylogenetic analysis revealed the presence of five distinct subpopulations within L. sativa with a variable degree of mixture across cultivar groups, confirming previous studies using both short-read genotyping-based techniques (Park et al., 2021; Park et al., 2022; Stoffel and van Leeuwen, 2012) and microsatellites (Rauscher and Simko, 2013). This could be related to the fact that in lettuce breeding, different horticultural types may be used in the pedigree scheme. In the collection assayed, we found a major clustering of butterhead genotypes when compared to the rest. This tendency contrasted with Park and colleagues (2021), who reported instead a greater separation of iceberg accessions from the other types in a collection of 441 individuals. Despite finding a slight differentiation according to geographical provenance, the effect due to the composition of the diversity panel assayed in terms of horticultural types and represented countries must be considered. Furthermore, for breeding and research materials, the reported origin often matches the places where the selection is carried out, thus providing an additional confounding effect. Several factors could affect the subpopulations enclosed in germplasm collections, such as the management practices occurring in the holding genebanks (e.g., level of heterozygosity retained and duplications), the areas of sampling of materials, or the biological status of accessions. Iceberg types investigated by Park et al. (2021) were mostly patented lines from the USDA, whereas we assayed mostly breeding materials, thus suggesting the presence of accessions still under development. The possibility to discover de novo polymorphisms free from any sequencing ascertainment bias and at affordable costs commensurable to other next-generation genotyping methodologies designates SPET as an efficient tool for population genomic analysis in lettuce.

The advances in genomics and cutting-edge genotyping technology have contributed to the growing availability of large-scale genotypic data of germplasm resources for various crops. The analysis of the genetic underpinnings of complex traits used in GWAS has benefited greatly from the ability to link phenotypic data to genomic sequence data. GWAS has proven to be an effective method for finding genetic variations that are significantly more common for a specific phenotype in unrelated individuals (Xiao et al., 2022). Owing to the greater number of recombination events occurring in natural populations, the advantage over bi-parental mapping populations depends on a larger genetic base to exploit and on higher map resolution (Han et al., 2020). Over the past years, the GWAS computing efficiency has been improved by developing different multivariate models that consider the family kinship inference and population structure covariates to enhance the power of associations and decrease the rate of false positives (Wang and Zhang, 2021). GWAS has been performed with the aim of investigating the potentiality of the SPET panel for candidate gene detection. To that end, we focused on four main agronomic traits driving the selection of cultivated lettuce cultivars and underlying market and consumer preferences. To test the most likely candidate regions underpinning the variation of the considered traits, different models were implemented. As expected, the GLM detected the highest number of STA in all traits, although this model accumulates several false positives, which are eliminated by incorporating additional correcting factors involving a multi-dimensional genome scan able to simultaneously estimate all marker effects (Wang et al., 2014: Chaurasia et al., 2021). By combining multivariate models, we identified seven candidate regions across chromosomes 5, 7, and 9 for the assayed traits. For seed color, the STA found on chromosome 7 confirmed a previous investigation reporting three associations in a 12-Mbp region spanning 69.87 Mbp to 80.63 Mbp (Kwon et al., 2013). We better refined the position at 50.40 Mbp near deah5 , an ATP-dependent RNA helicase involved in abscisic acid and stress responses in the acquisition of embryogenic competence (Almeida et al., 2020). The cytokinin dehydrogenase 3 detected within the association may regulate cell division as well as a large number of developmental events in plants (Schmülling et al., 2003). The two candidates may therefore play a role in seed coat development and color.

The position of STA located on chromosomes 5 and 9 for leaf color traits agreed with previous studies (Zhang et al., 2017; Su et al., 2020; Wei et al., 2021). On chromosome 5, Zhang et al. (2017) reported the lead SNPs for leaf color at less than 150 bp (86,123,627, 86,123,633, and 86,123,651) from the top association for leaf anthocyanin color. In the same study, the association on chromosome 9 was in the same region at 17.45 kb (152,892,248) from the top- ranked STA found in this study. These regions are reported to harbor two genes RLL2 (Red Lettuce Leaves 2) and ANS (Anthocyanin Synthase) that encode key enzymes for anthocyanin biosynthesis. We found four main candidate genes. Phototropin 2 and mlo like protein 11 both play a key role in leaf development and physiology. Phototropin 2 is primarily involved in the reception of light direction in the blade and has been demonstrated to promote leaf expansion and flattening (Legris et al., 2021). In Pistacia chinensis, Phototropin 2 has been reported to be involved in the signal transduction for anthocyanin accumulation during leaf coloration in autumn (Song et al., 2021), whereas in octaploid strawberry, it was involved in anthocyanin accumulation in strawberry fruits (Kadomura-Ishikawa et al., 2013).

The mlo like protein 11 is part of the large family of proteins that regulates pathogen defense and leaf cell death (Pozharskiy et al., 2022). No previous report indicates any function of mlo like protein 11 in leaf color. A general transcription factor ( 3C polypeptide 6 ) was found to be involved in outer leaf color on chromosome 9. In plants, transcription factors regulate secondary metabolism (Vom Endt et al., 2002) and are potential candidates for plant organ pigmentation (Ban et al., 2007; Zhou et al., 2014; Su et al., 2020).

The strongest signals found for bolting time at 164.43 Mb on chromosome 7 confirmed previous evidence. Indeed, several studies consistently supported the importance of chromosome 7 for lettuce flowering control (Kwon et al., 2013; Sthapit Kandel et al., 2020; Lee et al., 2021; Rosental et al., 2021). Despite the exact comparisons of the candidate region not always being possible, owing to the different marker system used (Han et al., 2021), our study supports whole-genome resequencing data findings (Wei et al., 2021), which detected a strong association at 164.5 Mbp in correspondence to phytochrome c involved in delaying of flowering. With the same effect, the strong STA found in the present study was near far1 (far-red impaired response 1), a component of the phytochrome A and putatively involved in regulating light control during the developmental stage (Siddiqui et al., 2016; Liu et al., 2020). far1 directly activates the expression of the evening gene ELF4 that plays a key role in the circadian flowering clock. In Arabidopsis, it negatively regulates flowering time in synergy with other FRS (FAR-Related Sequence) and FRF (FRS-Related Factor) genes (Ma and Li, 2018). The variation of far1 expression has also been reported to regulate shoot growth and flowering time in roses.

The creation of a novel SPET assay in lettuce was described in this work, and its potential for genetic diversity and GWAS research was demonstrated by comparing the results with earlier discoveries using different genotyping technologies. Additional research could weigh the benefits and drawbacks of SPET in comparison to whole-genome short and long read sequencing.