Genome‐wide association studies in a diverse strawberry collection unveil loci controlling agronomic and fruit quality traits

Abstract

Strawberries (Fragaria sp.) are cherished for their organoleptic properties and nutritional value. However, breeding new cultivars involves the simultaneous selection of many agronomic and fruit quality traits, including fruit firmness and extended postharvest life. The strawberry germplasm collection here studied exhibited extensive phenotypic variation in 26 agronomic and fruit quality traits across three consecutive seasons. Phenotypic correlations and principal component analysis revealed relationships among traits and accessions, emphasizing the impact of plant breeding on fruit weight and firmness to the detriment of sugar or vitamin C content. Genetic diversity analysis on 124 accessions using 44,408 markers denoted a population structure divided into six subpopulations still retaining considerable diversity. Genome‐wide association studies for the 26 traits unveiled 121 significant marker‐trait associations distributed across 95 quantitative trait loci (QTLs). Multiple associations were detected for fruit firmness, a key breeding target, including a prominent locus on chromosome 6A. The candidate gene FaPG1, controlling fruit softening and postharvest shelf life, was identified within this QTL region. Differential expression of FaPG1 confirmed its role as the primary contributor to natural variation in fruit firmness. A kompetitive allele‐specific PCR assay based on the single nucleotide polymorphism (SNP) AX‐184242253, associated with the 6A QTL, predicts a substantial increase in fruit firmness, validating its utility for marker‐assisted selection. In essence, this comprehensive study provides insights into the phenotypic and genetic landscape of the strawberry collection and lays a robust foundation for propelling the development of superior strawberry cultivars through precision breeding.

Abbreviations

• AsA
• ascorbic acid
• BLINK
• Bayesian-information and linkage-disequilibrium iteratively nested keyway
• CR
• call rate
• FarmCPU
• fixed and random model circulating probability unification
• GLM
• general linear model
• GWAS
• genome-wide association study
• IBS
• identity by state
• KASP
• kompetitive allele-specific PCR
• MLM
• mixed linear model
• NJ
• neighbor-joining
• PCA
• principal component analysis
• PG
• polygalacturonase
• QTL
• quantitative trait locus
• SNP
• single nucleotide polymorphism

INTRODUCTION

Strawberry (Fragaria × ananassa Duch.) fruit is highly consumed and appreciated for its vibrant red color, flavor, aroma, and nutritional value. These sensory attributes and potential health benefits to consumers have contributed to this fruit crop's economic importance, reaching a world production of 9.5 million tonnes in 2022 (FAOSTAT database, 2024; ). A substantial trait diversity is observed in strawberry varieties due to their elevated heterozygosity, a relatively brief breeding history spanning less than three centuries, and breeding for diverse climates and production systems (Darrow, 1966; Hardigan, Lorant et al., 2021). Breeding efforts in strawberry have been primarily directed toward preserving high yields, enhancing pest resistance, as well as improving fruit weight, firmness, and flavor (Senger et al., 2022; Verma et al., 2017; Whitaker et al., 2017). Among fruit quality traits, increasing fruit firmness has been of maximum interest, as it influences the organoleptic quality but, most importantly, the fruit's postharvest life and market value. During strawberry ripening, fruit firmness decreases dramatically due to the induction of enzymes involved in the degradation or modification of the cell wall, including polygalacturonases (PGs), pectin methylesterases, pectate lyases, endoglucanases, and expansins (reviewed in Moya-Leon et al., 2019). Extending postharvest life is essential for long-distance transportation. For instance, according to the FAOSTAT database (2024, ), Spain produced 325,880 metric tons of strawberries in 2022, with more than 85% of this production transported to Germany and other European countries via road transportation.

In response to consumer preferences, there is a growing emphasis on enhancing fruit organoleptic and nutritional quality in current breeding programs. Therefore, many strawberry breeders worldwide address the improvement of sugar/acid balance, aroma, fruit color, and antioxidant capacity or ascorbic acid (AsA) content (Mezzetti et al., 2016; Senger et al., 2022; Verma et al., 2017). On the other hand, the main breeding objectives for growers are focused on agronomic traits related to yield (number of fruits and size), resistance to pests and diseases (Phytophthora crown rot, Fusarium wilt, Anthracnose, or powdery mildew, among many others), and key physiological traits related to plant vegetative propagation (runnering time and number of runners), flowering habit, and fruit precocity.

The cultivated strawberry (Fragaria × ananassa Duchesne ex Rozier) is a complex allo-octoploid species (2n = 8× = 56) originating from the hybridization of two octoploid wild species: The South American Fragaria chiloensis (L.) and the North American Fragaria virginiana Miller (Darrow, 1966; Hancock, 1999). Their octoploid genome is formed by four subgenomes derived from four diploid progenitor species. Two of these subgenomes can be traced back to the ancestors Fragaria vesca L. and Fragaria iinumae Makino, whereas the donors of the remaining two subgenomes might be extinct and have not been clearly identified (Edger et al., 2018, 2019; Hardigan, Lorant et al., 2021; Liston et al., 2020; Tennessen et al., 2014). Genome-wide association studies (GWASs) were initially implemented in cultivated strawberry using the reference genome of the diploid F. vesca. Those first analyses led to the successful identification of genomic regions linked to several agronomic and fruit quality traits, such as disease resistance (Pincot et al., 2018), flowering time (Hardigan et al., 2018), or volatile organic compounds (VOCs) influencing fruit aroma (Barbey et al., 2021). Since the release of the first octoploid strawberry genome sequenced with chromosome-level resolution (Edger et al., 2019), association studies in this crop have benefited from subgenome-specific detection. These studies have explored the genetic control of fruit quality traits such as color (Castillejo et al., 2020), firmness (Hardigan, Lorant et al., 2021), or VOCs (Fan et al., 2022), among others. Despite all these studies, the number of traits explored by GWAS is still limited for this crop. It is worth noting that GWAS results could be influenced by various factors, including association models, population size and structure, phenotype, and environmental factors (Gai et al., 2023).

This research aims to uncover new loci contributing to the variation of important breeding targets and develop biotechnological tools for strawberry molecular breeding. To achieve this, we conducted GWAS using a population comprising, depending on the season, 95–124 different strawberry accessions from the Fragaria germplasm collection at Instituto Andaluz de Investigación y Formación Agraria y Pesquera (IFAPA, ESP138). This population included varieties of F. × ananassa and some hybrids with F. chiloensis, covering a broad geographic and temporal diversity, as well as a comprehensive representation of accessions adapted to the Californian-Mediterranean climate. Accessions were genotyped using the 50K FanaSNP Axiom array (Hardigan et al., 2020) mapped on the reference octoploid genome (Edger et al., 2019). As an initial step preceding the GWAS, we analyzed the similarity among accessions, the population's structure, and their genetic diversity. Over three consecutive seasons, we phenotyped this collection for 26 agronomic and fruit quality traits in order to evaluate the influence of environmental factors and identify loci that exhibit stability over the years. We report previously detected and novel quantitative trait loci (QTLs) for important breeding targets and other characters used as descriptors in the examination of distinctness, uniformity, and stability of the International Union for the Protection of New Varieties (UPOV). Among all those traits, we focused our attention on a major QTL for fruit firmness on chromosome 6A. We searched the haploblock in linkage disequilibrium (LD) for candidate genes and based on functional annotation, expression analyses, and previously reported FaPG1 studies (Garcia-Gago et al., 2009; Quesada et al., 2009), we propose FaPG1 as the underlying gene. Furthermore, we developed and validated a kompetitive allele-specific PCR (KASP) assay to assist in the selection of strawberries with higher fruit firmness.

Core Ideas

• A collection of 124 diverse strawberry accessions was phenotyped for 26 agronomic and fruit quality traits.
• Several quantitative trait locus (QTL) controlling agronomic and fruit quality traits were detected by genome-wide association studies.
• Natural variation in FaPG1 expression is associated with a major and stable QTL for fruit firmness.
• A marker assay was developed and validated for marker-assisted improvement of fruit firmness in strawberry.

MATERIALS AND METHODS

Plant material

The study was conducted during three consecutive seasons (2018–2019, 2019–2020, and 2020–2021). The initial experimental population consisted of 134 strawberry accessions from the Fragaria germplasm collection at IFAPA, including cultivated strawberry varieties (F. × ananassa) as well as hybrids with F. chiloensis to cover a wide diversity (Table S1). The cultivar Mieze Schindler, for which two different accessions from different sources (ACC_248 and ACC_307) were available in the collection, was used as an internal control. Young plants of each genotype were obtained each season after vegetative multiplication of accessions from the IFAPA Fragaria germplasm collection. Six rooted plants per genotype were transplanted during November to a shaded greenhouse in Málaga (36°40′22.5″ N 4°30′13.2″ W), Spain, following a randomized three block design. To simulate field conditions, the experimental population was cultivated in 1 m × 1 m containers filled with a mixture of universal substrate with river sand in a 3:1 v/v ratio. Eight plants (four accessions × two replicates) were planted in each container at a distance of 25 cm, which is similar to the strawberry planting density in commercial fields at Huelva (37°14′26.8″ N 6°48′06.3″ W), Spain. Plants were irrigated automatically using a drip system 3 days per week, increasing to 7 days during summer. During all seasons, plants were treated monthly against white spider mite and red spider mite, cleaned weekly of old leaves and weeds, and fertilized with different concentrations of nitrogen (20–30 units/ha), phosphorus (7–12 units/ha), potassium (11–44 units/ha), calcium (6–10 units/ha), magnesium (2–4 units/ha), and organic matter (10–30 units/ha) according to the time of year. Plants were regularly treated against powdery mildew during the second and third seasons.

DNA extraction, genotyping, and pedigree confirmation

Genomic DNA was extracted using 40–50 mg of freeze-dried young leaf tissue following the instructions of the Omega DNA extraction kit (Omega Bio-tek; ) with minor modifications specified in Hardigan et al. (2018) and quantified in a Qubit fluorometer (Thermo Fisher Scientific). DNAs from individuals were submitted to Segalab for genotyping using the 50K FanaSNP Axiom array (Hardigan et al., 2020). The generated data files were analyzed with the Axiom Analysis Suite v.5.1.1.1 software (Thermo Fisher Scientific) following the Best Practices Genotyping Workflow described in the software documentation. Markers were filtered for a call rate (CR) ≥97%, and the parameters established for polyploid species by default. Single nucleotide polymorphism (SNP) and indel markers assigned into the “poly high resolution” (polymorphic with three defined clusters) or “no minor homozygote” (polymorphic with two defined clusters) classifications were selected for further analyses. Samples with a CR >97% were obtained, although two accessions with CRs of 93% and 94% were manually selected and retained.

Similarity between accessions was analyzed using the identity by state (IBS) distance matrix using PLINK v.1.9 software (Chang et al., 2015). Accessions with higher similarity than that obtained between the two accessions of Mieze Schindler (ACC_248 and ACC_307; IBS > 99.78%) were considered propagation errors or synonyms. In order to determine one or the other situation, all genotypes were compared with data obtained for those same accessions in a previous genotyping of the IFAPA collection performed with an independent DNA extraction and the Affymetrix IStraw90K Axiom array (Bassil et al., 2015). A total of 10 accessions that turned out to be propagation errors, duplicates and/or synonyms were removed from further analyses, remaining a total of 95, 114, and 124 accessions in the first, second, and third seasons.

Genetic diversity and population structure

Genetic diversity analyses were performed with the 44,408 high-quality SNPs and indels selected from the 50K FanaSNP Axiom array (Hardigan et al., 2020) in 124 accessions of the experimental population using different methods. Principal component analysis (PCA) was implemented using TASSEL v.5.0 (Bradbury et al., 2007). The population structure of the experimental population was analyzed using 44,002 SNPs (w/o indels) and the software STRUCTURE v.2.3.4 (Pritchard et al., 2000). Two to 20 subpopulations were evaluated with an admixture model, 10,000 burn-in steps, 100,000 Markov-Chain Monte Carlo (MCMC) steps, and 10 replicates per K-values. This analysis was carried out at the facilities of the Supercomputing and Bioinnovation Center of the University of Malaga () using StrAuto (Chhatre & Emerson, 2017). The optimal number of subpopulations (K-value) was established using STRUCTURE HARVESTER (Earl & vonHoldt, 2012) and the Evanno method (Evanno et al., 2005). Sample orders were calculated using CLUMPP v. 1.1.2 (Jakobsson & Rosenberg, 2007), and cluster results were visualized using Distruct v.1.1 (Rosenberg, 2004). Neighbor-joining (NJ) phylogenetic tree was constructed with TASSEL v.5.0 (Bradbury et al., 2007) from the 1-IBS distance matrix. The dendrogram obtained was drawn and visualized with Fig Tree v.1.4.4 (). Finally, the pairwise fixation index (F_ST) was calculated using the Hierfstat v0.5-7 tool for R (Goudet, 2005) with data obtained from the STRUCTURE, to know the degree of population structuring.

Phenotypic evaluation

Phenotypic characterization for a total of 26 morphological, physiological, and fruit quality traits was carried out during three consecutive years. Number of evaluated traits varied from 15 in the first season to 25 in the second and third seasons (Table 1). Trait data collection was carried out in the greenhouse weekly during the optimum of plant and/or fruit development. Traits were noted by the letters FL if they relate to flowers, LE if they relate to leaves, PL for plant-related, RU for runners, and FR for fruit-related, followed by the name of the traits or by an abbreviation of it. Reference scales established for variety protection by UPOV ( Table 1) were used for 12 traits, which were evaluated in the six replicates from each accession.

TABLE 1 Summary of the traits evaluated during the three seasons.

Trait	2020–2021	2019–2020	2018–2019	Correlations
Name	Abbreviation	Units	Mean ± SD	Range	Mean ± SD	Range	Mean ± SD	Range	2020–2021 vs. 2019–2020	2020–2021 vs. 2018–2019	2019–2020 vs. 2018–2019
Flowering time	FL_Days	day	111.95 ± 36.11	48–220	164.62 ± 38.63	58–237	111.56 ± 34.32	27–200	0.49*	0.34*	0.33*
Green color in petals	FL_GCP	2-point scale	1.14 ± 0.35	1– 2	1.02 ± 0.06	1–2	1.03 ± 0.12	1–2	0.14	0.15	0.15
Red color in petals	FL_RCP	3-point scale	1.13 ± 0.34	1–3	1.03 ± 0.12	1–3	1.15 ± 0.37	1–3	0.50*	0.48*	0.41*
Flower diameter	FL_Diameter	3-point scale	5.00 ± 0.32	3–7	4.95 ± 0.60	3–7	4.80 ± 0.83	3–7	0.23*	0.13	0.44*
Leaf color upper side	LE_CUS	5-point scale	3.27 ± 0.68	1–5	3.17 ± 0.62	1–5	2.98 ± 0.33	1–5	0.57*	0.47*	0.35*
Leaf glossiness	LE_Glossiness	3-point scale	2.00 ± 0.47	1–3	2.13 ± 0.45	1–3	1.97 ± 0.35	1–3	0.40*	0.41*	0.49*
Leaf size	LE_Size	3-point scale	4.74 ± 0.99	3–7	4.94 ± 0.95	3–7	–	–	0.59*	–	–
Leaf mildew	LE_Mildew	5-point scale	–	–	–	–	0.34 ± 0.76	0–4	–	–	–
Plant height	PL_Height	cm	19.01 ± 3.98	7.17–31.50	16.09 ± 3.77	6.40–27.83	14.38 ± 3.86	4–23.50	0.63*	0.61*	0.60*
Plant diameter	PL_Diameter	cm	33.22 ± 8.35	11.50–60.67	28.60 ± 6.70	14.40–47.50	25.42 ± 5.63	14–38.00	0.65*	0.55*	0.59*
Runnering time	RU_Days	day	223.99 ± 22.77	121–264	247.89 ± 22.81	202–303	256.17 ± 29.68	175–303	0.67*	0.82*	0.71*
Runner number	RU_Num	No. of runners	5.83 ± 4.33	0–20	4.88 ± 3.60	0–19	3.14 ± 2.68	0–12	0.72*	0.66*	0.71*
Runner anthocyanins	RU_Anthocyanins	5-point scale	2.79 ± 0.87	1– 5	2.33 ± 0.81	1–5	2.62 ± 0.97	1–5	0.53*	0.78*	0.73*
Fruit number	FR_Num	No. of fruits	15.14 ± 6.10	4–32	9.03 ± 4.54	0–25	6.49 ± 3.36	0–14	0.60*	0.22*	0.15
Fruit weight	FR_Weight	g	9.33 ± 5.29	0.67–23.44	8.04 ± 4.55	0.83–27.11	–	–	0.71*	–	–
Fruit width of band without achenes	FR_WB	5-point scale	3.28 ± 1.07	1–9	3.12 ± 1.01	1–9	2.84 ± 0.97	1–9	0.53*	0.61*	0.59*
Fruit external color	FR_EC	7-point scale	4.67 ± 0.93	1–7	4.76 ± 0.80	1–7	–	–	0.64*	–	–
Fruit external color:Lightness	FR_L*	black (0)–white (100)	30.18 ± 3.89	23.30–50.19	30.95 ± 3.86	23.23–52.38	–	–	0.77*	–	–
Fruit external color:Red-green	FR_a*	green (−a)–red (+a)	38.60 ± 4.40	9.85–49.40	41.27 ± 4.39	9.59–47.85	–	–	0.78*	–	–
Fruit external color:Yellow-blue	FR_b*	blue (−b)–yellow (+b)	20.89 ± 3.52	13.14–29.29	22.78 ± 3.66	14.63–32.29	–	–	0.57*	–	–
Fruit internal color:Flesh	FR_ICF	6-point scale	4.59 ± 0.98	1–6	4.14 ± 1.31	1–6	–	–	0.44*	–	–
Fruit internal color:Core	FR_ICC	3-point scale	2.09 ± 0.54	1–3	2.01 ± 0.66	1–3	–	–	0.36*	–	–
Fruit firmness	FR_Firmness	g	170.91 ± 54.21	100–303.13	197.27 ± 67.24	101.25–340	211.36 ± 50.48	127.5–324.38	0.87*	0.27*	0.20*
Fruit acidity	FR_Acidity	g kg−1	11.80 ± 2.50	4.70–19.50	11.30 ± 2.40	1.20–20.10	–	–	0.68*	–	–
Fruit soluble solids content (SSC)	FR_Brix	%	7.61 ± 1.58	4.93–11.80	8.09 ± 1.13	5.87–10.90	–	–	0.65*	–	–
Fruit ascorbic acid	FR_AsA	mg kg−1	473.2 ± 116.6	161.1–834.9	478.6 ± 126.3	208.7–812.5	–	–	0.73*	–	–

Powdery mildew natural infection on the plant (no inoculation performed) was evaluated only during the first season from November until August. Symptoms were evaluated weekly in the six replicates per accession using a 1–5 rating scale based on the presence of powdery mildew symptoms on the leaf: (0) healthy plant with 0% presence, (1) <10%, (2) 11%–25%, (3) 26%–50%, or (4) >50% leaf area affected.

Besides using the UPOV scale, external fruit color was measured in a total of eight fruits per accession with a chromameter (Chroma Meter CR 410, Konica-Minolta) using CIELAB color space values a* (green–red spectrum), b* (blue–yellow spectrum), and L* (brightness–darkness). Fruit weight was evaluated as the mean weight of 10 fully ripe fruits at the peak of production. Flowering and runnering time was recorded as the average date of the first open flower or runner in the six replicates of each accession, respectively.

Fruit firmness was evaluated using a penetrometer (Facchini Fruit Pressure tester FT 02) with a 3.5-mm probe in eight fruits per genotype, taking two measurements on each fruit. Measurements were taken at the same time of the morning during the production peak of each accession, which took place from April 1 to July 31 depending on the accession.

For quality traits evaluated in processed fruit, at least 25 commercially mature fruits of each genotype were harvested weekly, at approximately the same time of the morning, and during the production peak of the season. These fruits were immediately cut, frozen in liquid nitrogen, and stored at −80°C. Fruits of each accession were divided into three biological replicates consisting of a minimum of eight fruits. The set of fruits from each replicate was powdered in liquid nitrogen using a coffee grinder and stored at −80°C until analysis. The first fruits produced by the plant were discarded because they were not representative of each variety, as well as misshaped fruits. Fruit acidity (citric acid equivalents, g kg⁻¹) expressed on a fresh weight basis was measured with an automatic titrator (TitroLine easy, Schott Instruments, GmbH) by diluting 1 g of fruit powder from each replicate in 100 mL of distilled water. Titrations were carried out to a final pH of 8.1 using 0.01 N NaOH. Soluble solids content (SSC, %) was quantified using a digital refractometer (ATAGO Brix digital refractometer PR-32α) and a drop of blended fruit from the three replicates. AsA content (mg kg⁻¹) expressed on a fresh weight basis was quantified spectrophotometrically using the protocol of Fenech et al. (2021) with modifications: approximately 100 mg of the fruit powder from each replicate was homogenized in 1 mL of cold extraction buffer (3% metaphosphoric acid, 1 mM EDTA) and centrifuged at 14,000 rpm for 20 min at 4°C. The supernatant was filtered through a 0.45-µm nylon membrane. Absorbance measurements were carried out at 265 nm in the plate reader on UV-transparent 96-well plates (Greiner UV-Star 96-Well Microplate) after 30 min incubation of 20 µL of the filtered sample with 100 µL of phosphate buffer (0.2 MKH₂PO₄). After the measurement, 5 µL of the ascorbate oxidase enzyme (40 U mL⁻¹ in phosphate buffer) was added to each well, mixed, and incubated at room temperature for 20 min, and the absorbance was measured again. The AsA content was calculated by interpolating the difference of the absorbance on a standard curve with known concentrations of AsA (Sigma).

Correlations for each trait between the 3 years and correlations between the log-transformed data (log₁₀) of the 25 traits measured in the third season (2020–2021) were calculated using Pearson correlation in R version 4.3.2 using “stats” and “corrplot” packages. PCA of 124 accessions of the third season was performed with the log-transformed data (log₁₀) of the 25 traits in R using “FactoMineR” and “Factoextra” packages. Boxplots of fruit firmness across subpopulations were calculated using the R package “ggpubr.”

Genome-wide association study

GWAS was performed with the genotypic data from the 50K FanaSNP Axiom array (Hardigan et al., 2020) and phenotypic data from each season separately. Missing genotype data were imputed with Beagle v.4.0 (Browning & Browning, 2007), obtaining a total of 124 accessions genotyped with 44,196 markers mapped to the reference ‘Camarosa’ v1.0. a2 genome (T. Liu et al., 2021). For phenotypic data, the mean of the biological replicates for each accession in each year was used. Those traits that deviated from normality were transformed using the Cox–Box transformation (Box & Cox, 1964) only when transformation improved the normality. GWAS was conducted in the R package GAPIT v.3 (Lipka et al., 2012; Wang & Zhang, 2021) using the general linear model (GLM), mixed linear model (MLM), fixed and random model circulating probability unification (FarmCPU), and Bayesian-information and linkage-disequilibrium iteratively nested keyway (BLINK) analysis models (Huang et al., 2019; X. Liu et al., 2016; Price et al., 2006; Yu et al., 2006). All associations with p-value < 0.05 after the False Discovery Rate (FDR) control procedure were considered significant. Those significant markers that were within a distance of <5 Mb were considered a single QTL. To consider the population stratification, the genomic relationship matrix (G) was generated using the VanRaden algorithm (VanRaden, 2008) and included as random effect, while the top three principal components, explaining 26.1% of the variance, were used as covariates in the GWAS model. Haplotype blocks were computed with “SnpStats” package in R using ld function to calculate pairwise LD (D') and plotted with ld.plot function from the “gaston” package.

In silico candidate gene identification

Genes associated with the positions of significant SNPs were obtained by searching for the position of the significant SNP marker in the Genome Database for Rosaceae (GDR; ). The gene located at that specific position was annotated as candidate gene, or in cases where no gene was found, the flanking genes to it. Subsequently, a basic local alignment search tool analysis was performed to determine the function assigned to each of the associated genes using the ‘Camarosa’ reference genome v1.0 (Edger et al., 2019) with the new a.2 annotation (T. Liu et al., 2021).

Candidate genes within the haploblock for fruit firmness on chromosome 6A were manually searched based on functions related to cell-wall degradation. Previously, annotations of genes located within the QTL interval were retrieved from ‘Camarosa’ reference genome assembly v1.0 a.2 (Edger et al., 2019). Expression level of candidate genes during ‘Camarosa’ fruit ripening and in ripe ‘Senga Sengana’ fruit was obtained from RNAseq data from previous studies (Castillejo et al., 2020; Sánchez-Sevilla et al., 2017) but here mapped to either the ‘Camarosa’ v1.0 a.2 (Edger et al., 2019) or the ‘Royal Royce’ FaRR1 (Hardigan, Feldmann et al., 2021) genomes using the Tuxedo suit as previously described (Castillejo et al., 2020).

RNA isolation and gene expression analysis

We generated two pools of fruits selecting 30 accessions contrasting in fruit firmness among the 124 accessions used in the third-season (2020–2021) GWAS. Each pool in triplicate consists of an equivalent amount of fruit tissue stored at −80°C from each of the three biological replicates of 15 accessions. We isolated total RNA from the three biological replicates from each pool of 15 accessions (high fruit firmness versus low fruit firmness) using the Plant/Fungi Total RNA Purification Kit (Norgen Biotek). Prior to reverse transcription, RNA was treated with DNAse Turbo (Invitrogen) following the manufacturer's instructions. A total of 800 ng of RNA was used for retrotranscription using the High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems by Thermo Fisher Scientific). Primers for quantification of FaPG1 by qRT-PCR, FaPG1_F: GCTCCTGGTGACTTTGATGT and FaPG1_R: ACTCTACTTGGCGTTGTTGC, were described in López-Casado et al. (2023). Subgenome-specificity of primers for FaPG1 was confirmed by comparison to the other PG homoeologs and the most similar homologs retrieved from the GDR database (). qRT-PCRs were conducted in a CFX96 Real-Time System (Bio-Rad) using 2 µL of a 1:5 dilution of the cDNA, 250 nM of each specific primer, and the TB Green Premix Ex Taq (Takara) in a final volume of 10 µL. Three technical replicates were used for each biological replicate. The polymerase chain reaction (PCR) program involved an initial denaturation at 95°C for 30 s, followed by 39 cycles of denaturalization at 95°C for 10 s, and annealing at 60°C for 25 s. Relative expression of the FaPG1 gene was calculated by the 2^−ΔΔCT method, using FaGAPDH and FaDBP as reference genes (Muñoz-Avila et al., 2022) and the high-firmness pool as the reference sample. Statistical analysis was performed using GraphPad Prism 8.0 software. The normality of data distribution was evaluated by Kolmogorov–Smirnov and then analyzed by t-test.

Kompetitive allele-specific PCR assay

A KASP assay for AX-184242253, one of the significant SNPs detected on chromosome 6A, was designed using PolyOligo () and the ‘Camarosa’ v.1.0 reference genome. To further reduce the possibility of off-targets, primers were manually adjusted for increased sub-genome specificity (FIRM6A-KASP_REF: GAAGGTCGGAGTCAACGGATTattgccacaaatgcaagagG, FIRM6A-KASP_ALT: GAAGGTGACCAAGTTCATGCTattgccacaaatgcaagagA, and FIRM6A-KASP_COM: GCTAAAACATGACAACATCATCTTAC). KASP assays (Semagn et al., 2014) were performed on a CFX96 real-time thermal cycler (Bio-Rad) using the KASP-TF Master Mix (LGC Genomics). The PCR program consisted of an activation at 94°C for 15 min, 10 touchdown cycles of 94°C for 20 s and 61°C–56°C for 1 min (decreasing 0.5°C per cycle), followed by 28 cycles of 94°C for 20 s and 56°C for 1 min, and plate reading after one step of 1 min at 37°C. One or two recycling programs consisting of three additional cycles were necessary for cluster improvement and/or amplification of some accessions. Statistical analyses were performed using Student's t-test implemented in GraphPad Prism version 8.0 software.

RESULTS

Genetic diversity and population structure

After the analysis of similarity using the IBS calculations and pedigree confirmation of the 134 accessions of the experimental population, 10 of them were removed from the study due to possible propagation errors (five accessions), duplicates and/or synonyms (Table S2). Among synonyms, accessions Revada, Georg Soltwedel, Africa, and Gorella White Pulp were considered synonyms and only Africa (ACC_256) was retained in the studies. Similarly, Victorian Nameless and Surprise des Halles were considered synonyms and only Surprise des Halles (ACC_266) was retained. Therefore, genetic diversity, population structure, and F_ST index were analyzed using a total of 124 accessions and 44,408 high-quality SNPs and indels.

Evaluation of population structure using PCA resulted in the two first principal components (PCs) explaining 16% and 5% of the variance (Figure 1A). PCA organized the experimental population by time and geographic adaptation, with the most modern Californian-Mediterranean varieties on the left side (in red color) and another distinct group including the oldest European varieties adapted to northern territories (in blue color) on the right side. The rest of varieties were distributed between these two groups according to time and climatic/geographical adaptation. Recently developed hybrids between Californian varieties and F. chiloensis (highlighted in green) clustered in the lower and central parts of the graph.

[IMAGE OMITTED. SEE PDF]

To assess the level of genetic stratification, we used the software STRUCTURE, which showed a diversification of the strawberry collection into six subpopulations (Figure 1B,C). The largest subpopulation includes classic Californian cultivars ‘Chandler’ and ‘Camarosa’ and more recent Californian/Mediterranean/Floridian varieties, such as ‘Candonga’ or ‘Fortuna’ (subpopulation 1; Figure 1C). The second subpopulation (in dark purple in Figure 1C) contains a considerable level of admixture and is formed mainly by older Californian/Floridian varieties, such as ‘Pajaro’, ‘Toro’, or ‘Aiko’, and some Japanese accessions such as ‘Aiberry’ or ‘Toyonoka’. The oldest European accessions from northern territories belong to a third subpopulation (in blue in Figure 1C), which includes ‘Laxton’s Noble’ and ‘Vicomtesse’. Two subpopulations act as a transition between the blue and purple (subpopulations 4 and 5; Figure 1C). In subpopulation 4 (light purple), stand out some varieties defined by Darrow (1966) as “large-fruited strawberries best adapted to tropical climates” (‘Blakemore’, ‘Florida Ninety’, and ‘Missionary’) and classic North-American varieties as Howard 17. Subpopulation 5 (in pink) contains only two accessions from the Netherlands, ‘Tenira’ and ‘Gorella’. Finally, the population structure analysis discriminated a sixth subpopulation consisting of varieties introgressed with F. chiloensis (in green; Figure 1C) that cluster with other related hybrids belonging to subpopulation 1.

Despite the identification of six subpopulations, the degree of stratification was scarce, and subpopulations overlapped with each other in the dendrogram. This was corroborated by a fixation index value (F_ST) of 0.1343 for the experimental population, indicating a limited differentiation between subpopulations and the absence of clearly defined clusters in the collection.

Phenotypic variation of agronomic and fruit quality traits within the strawberry germplasm collection

The phenotypic evaluation of the experimental population was carried out during three consecutive seasons for morphological, physiological, and fruit quality traits. In the three seasons, a large phenotypic variation was observed among accessions for all evaluated traits (Table 1). This variation was consistent among the 3 years for most of the traits, with a medium-high correlation coefficient (0.40–0.80) between seasons, except for the green color of petals, which was not significant between any season, suggesting that it is a trait with a strong environmental influence (Table 1). For flower diameter and fruit number, high correlations were only observed between the first and second or the second and third years, respectively. High correlations among the three seasons were observed for runnering time, number, and color, suggesting a much higher genetic control than the genotype × environment interaction (G × E) and the environmental effects on these traits. Fruit firmness ranged from 1 to 34 N during the three seasons and showed a 0.87 correlation between the second and third seasons, but a low correlation with the first season. In contrast, the width of fruit band without achenes displayed high correlations between the 3 years, again indicating a lower environmental influence.

To identify relationships between the traits, we analyzed Pearson correlations between them using data from the third season (2020–2021; Figure 2A; Table S3). A total of 43 positive and 35 negative significant correlations were observed, of which 11 were considered high (±0.60 to ±0.80), 14 were moderate (±0.40 to ±0.60) and most of the significant correlations were, therefore, low (<±0.40). In general, we observed medium-high positive correlations between traits related to plant vigor (LE_Size, PL_Height, and PL_Diameter) and fruit or runner production (FR_Num, FR_Weight, and RU_Num). Among them, the highest significant correlations were observed between leaf size with plant height and plant diameter (LE_Size/PL_Height: +0.67 and LE_Size/PL_Diamenter: +0.66), and between plant height and plant diameter (PL_Height/PL_Diameter: +0.78). On the other hand, runnering time (RU_Days) was negatively correlated with four traits: PL_Height, PL_Diameter, LE_Size, and number of runners (RU_Num), with the latter correlation standing out with a value of −0.72. Fruit traits related to color also showed medium and high correlations. We observed positive correlations between the two internal fruit color traits (FR_ICF/FR_ICC: +0.65) and between fruit brightness and yellow hue (FR_L*/FR_b*: +0.77). Negative correlations have been observed between external color (visual scale) and brightness (FR_EC/FR_L*: −0.74) or yellow hue (FR_EC/FR_b*: −0.69). This is to be expected since higher values of L* or b* are associated with lighter fruit colors (Zorrilla-Fontanesi et al., 2011). A strong positive correlation (0.66) was observed between FR_Firmness and FR_Weight. Strikingly, FR_Firmness was negatively correlated with fruit acidity (FR_Acidity), FR_Brix, and FR_AsA content (−0.57, −0.55, and −0.38, respectively). A low positive correlation (0.37) was observed between FR_Brix and FR_Acidity. As previously reported, FR_Weight was negatively correlated with FR_Brix, FR_Acidity, and FR_AsA (−0.49, −0.48, and −0.34, respectively; Hummer et al., 2023; Muñoz et al., 2023; Stevens et al., 2007; Wada et al., 2020). Albeit weakly, fruit number was also negatively correlated with FR_Brix and FR_AsA.

[IMAGE OMITTED. SEE PDF]

A PCA of phenotype data supported the observed relationships (Figure 2B,C). For example, traits related to plant vigor cluster together in the loadings plot and had an important contribution to the separation of strawberry accessions across PC2, together with runner number and runnering time, the latter in the opposite direction (Figure 2B). ACC_1062, at the right top of the plot, is an example of the most vigorous genotypes (Figure 2C). Similarly, strawberry accessions were separated across PC1 based on fruit color, fruit firmness, and weight, which were also great contributors. Accessions 845, 872, and 751, on the left side, produce white or light-red fruits, while accessions with the darkest fruits appear on the opposite side of the PCA (Figure 2C). These three traits separated older from recent accessions across PC1 in accordance to their population structure. For instance, fruit firmness in the third season showed significant differences between the six subpopulations, indicating a strong selection over this trait across breeding history (Figure S1). Another important contributor to PC1, but in the opposite direction to fruit firmness, was FR_Brix, in accordance with the significant negative correlation observed between them.

Identification of phenotype-genotype associations for agronomic and fruit quality traits

Association analyses have been carried out for each year separately with the GLM and MLM single-locus test models and with the FarmCPU and BLINK multi-locus test models. A total of 121 significant marker-trait associations (FDR-adjusted p-value < 0.05) were detected in 19 of the 26 evaluated traits, distributed in 95 QTLs, in the three seasons and with the four models used (Table 2; Table S4). The number of associated SNPs varied largely between seasons and models (Table 3). In the first, second, and third seasons, nine, 36, and 105 marker/traits associations were detected, respectively. Regarding GWAS models, the majority of associations were detected with FarmCPU and BLINK (51 and 56), next 36 using GLM, and only seven with MLM. In general, SNPs associated with a trait detected with one model were also detected above the background with other models albeit not being significant.

TABLE 2 Number of significant SNPs and QTLs detected in the GWAS for the 19 traits in any of the three seasons.

Trait	Significant SNPs	Significant QTLs
PL_Diameter	2	2
RU_Anthocyanins	1	1
RU_Days	10	8
RU_Num	10	8
LE_CUS	4	4
LE_Glossiness	1	1
LE_Mildew	1	1
FL_RCP	10	8
FL_Diameter	1	1
FR_Num	12	12
FR_Weight	4	4
FR_WB	1	1
FR_L*	6	5
FR_a*	13	12
FR_b*	4	4
FR_EC	4	4
FR_Acidity	2	2
FR_Brix	10	10
FR_Firmness	25	7
Total	121	95

TABLE 3 Summary of significant SNPs detected in each season and with each of the four GWAS models.

Model	2020–2021	2019–2020	2018–2019	Total SNPs /model
GLM	29	6	1	36
MLM	6	1	0	7
FarmCPU	31	17	3	51
BLINK	39	12	5	56
Total SNPs/season	105	36	9	150

No QTL was detected for flowering time (FL_Days), petal greening (FL_GCP), leaf size (LE_Size), plant height (PL_Height), fruit internal color (FR_ICF and FR_ICC), or AsA (FR_AsA). Out of the total 95 QTLs, only six were detected in more than one season (Table S4). Only one QTL was detected in one individual season for the presence of anthocyanins in the runner (RU_Anthocyanins), for leaf glossiness (LE_Glossiness), leaf powdery mildew (LE_Mildew), flower diameter (FL_Diameter), and for presence of an achene-free band on the fruit (FR_WB; Figure 3; Figure S2). Two independent QTLs on chromosomes 6D and 7A were detected for fruit acidity (FR_Acidity; Figure 4A) and another two for plant diameter (PL_Diameter; Figure S3). Four QTLs were detected for external fruit color evaluated according to the UPOV visual scale (FR_EC; Figure 4B), the blue-yellow value of fruit external color (FR_b*; Figure S4), leaf color (LE_CUS; Figure S5), and for fruit weight (FR_Weight; Figure S6). A total of five QTLs were detected for the external fruit brightness (FR_L*; Figure S7). In contrast, seven or more QTLs were detected for the remaining seven traits (Table 2), indicating a more complex genetic architecture for those traits.

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

Runnering is an important trait for plant breeders and producers, and our results support a complex architecture, as eight QTLs were detected for both runnering time (RU_Days; Figure 5A) and number of runners (RU_Num; Figure 5B). One QTL for RU_Days on chromosome 2B and two for RU_Num on 4A and 5A were detected in two seasons (Figure 5; Table S4). Interestingly, the QTL on chromosome 5A has pleiotropic effects on both related traits, and the common SNP AX-184418817 delayed runnering by 18 days and decreased runner number by 1.6–3.3. In contrast to runnering, the presence of red coloration on petals (FL_RCP) is not currently a breeding target, and the genetic control of this trait has not been studied previously in strawberry. We have detected 10 associations in eight chromosomal regions considering all three seasons (Figure S8; Table S4), being the QTL on chromosome 3B common in the first and second seasons.

[IMAGE OMITTED. SEE PDF]

Crop production or yield is considered the most important breeding target and depends on both fruit number and weight. While only four QTLs were detected for FR_Weight during the third season, a total of 12 QTLs with small effects were detected for fruit number (FR_Num) in the 2019–2020 or the 2020–2021 seasons (Figure S9; Table S4). None of them was detected in two seasons. A total of 12 QTLs were also detected for the green-red value of fruit external color (FR_a*), with two of them, on chromosomes 1B and 6A, being detected in the last two seasons (Figure S10; Table S4). Remarkably, the QTL on 1B, at position 15.39 Mb, was shared with other related color traits (FR_L* and FR_EC), and one of the leading SNPs, AX-184967514 (Figure 4B; Tables S4 and S5), is located at the transcription factor FaMYB10, a positive regulator of anthocyanin biosynthesis and strawberry fruit color (Castillejo et al., 2020; Medina-Puche et al., 2014).

Our GWAS also showed a complex genetic architecture for SSC (FR_Brix), with one and nine significant associations in the second and third seasons, respectively (Figure S11; Table S4). Fruit firmness was the trait with the highest number of significant associations, with 25 SNPs distributed in a total of seven QTLs on chromosomes 1A, 2B, 3A, 5A, 5D, 6A, and 7A. Only one QTL was detected in the second season and the rest in the third (Figure 6; Table S4). QTLs on chromosomes 3A and 6A were each characterized for having 10 significant SNPs in an interval of 1.2 Mb and 4.3 Mb, respectively. Significant SNPs on 6A were detected with the four GWAS models. It is also worth noting that the effects on fruit firmness were about 0.3 N for each of these two QTLs (Table S4).

[IMAGE OMITTED. SEE PDF]

Linkage disequilibrium haploblock and candidate genes for fruit firmness QTL on chromosome 6A

As firmness is a vital breeding target, and due to the large effect of the QTL for FR_Firmness in chromosome 6A, we decided to focus our search for the underlying gene on this region. Besides, a QTL for firmness in the same region has been recently detected in other GWAS using wild and domesticated cultivars (Hardigan, Lorant et al., 2021) and a multiparent population (Cockerton et al., 2021), suggesting that it is conserved in diverse genotypes and environments. The majority of significant SNPs detected on 6A in our study spanned a chromosomal region of 0.89 Mb in high LD flanked by markers AX-166518257 and AX-184933268 (Table S4; Figure S12). This region contains 177 annotated genes in the ‘Camarosa’ reference genome, including four closely located genes with similarity to PGs (PGs: FxaC_21g15450, FxaC_21g15750, FxaC_21g15770, and FxaC_21g15780) and one transcript with similarity to endoglucanases (FxaC_21g15730; Table S6). Fruit softening during ripening is caused by disassembly of the cell wall and the dissolution of the middle lamella, and a set of enzymes, including PGs, mediates this process (Moya-Leon et al., 2019). Therefore, these genes are good candidates for controlling the observed natural variation on fruit firmness in the GWAS population.

We next studied the expression of these five genes during ‘Camarosa’ fruit ripening using RNAseq data from a previous study (Sánchez-Sevilla et al., 2017) but here mapped to either the ‘Camarosa’ v1.0 a.2 (Edger et al., 2019) or the ‘Royal Royce’ FaRR1 (Hardigan, Feldmann et al., 2021) genomes. Only one of them (FxaC_21g15770) showed significant expression in the analyzed tissues (Table S7). We also studied the expression levels of the candidate genes in ripe fruits of ‘Senga Sengana’, a cultivar with much softer fruits than ‘Camarosa’, using RNAseq data previously reported (Castillejo et al., 2020). Similar to ‘Camarosa’ fruits, FxaC_21g15770 was the only gene with significant expression (Table S7). Interestingly, this transcript encodes for the functionally characterized FaPG1 gene (AF380299; Quesada et al., 2009). As in previous reports, in our RNAseq analysis, FaPG1 was induced during fruit ripening, with the highest expression in the turning and ripe stages of the receptacle tissue. Although not being fully comparable due to being independent experiments, FaPG1 expression was 40–60 times higher in ripe fruits of ‘Senga Sengana’ than in ‘Camarosa’. While the gene annotation of FaPG1 in ‘Camarosa’ (FxaC_21g15770) contains the four expected exons, in the FaRR1 reference genome its annotation (Fxa6Ag103973) is not correct. The Fxa6Ag103973 gene spans eight predicted exons from two incomplete tandemly arranged PG genes. However, the two complete open reading frames of these two genes, FaPG1 and the next PG (FxaC_21g15750 in ‘Camarosa’), could be manually found in the ‘Royal Royce’ reference genome. An alignment of the four PG sequences from ‘Camarosa’ and ‘Royal Royce’ is shown in Figure S13. The gene FaPG1 was identical in ‘Camarosa’ and ‘Royal Royce’. Similarly, high similarity was observed between the ‘Camarosa’ and ‘Royal Royce’ alleles for the other three PGs, which were not significantly expressed in the analyzed tissues. FaPG1 nucleotide sequences from ‘Camarosa’ or ‘Royal Royce’ displayed six SNPs with AF380299, the FaPG1 gene from ‘Chandler’ (Quesada et al., 2009), and only three amino acid (Aa) differences in the protein sequence. None of the three Aas were located in any of the four conserved domains associated with PGs (Chen et al., 2016; Figure S14), suggesting that they encode active enzymes.

As FaPG1 was the only candidate gene being expressed, and to confirm its role in natural variation in fruit firmness, we quantified its expression levels by qRT-PCR in ripe fruit from two pools of accessions contrasting on fruit firmness (Figure 7). Expression of FaPG1 was 13-fold higher in the accessions with lower firmness, indicating that this is the most plausible cause of the observed variation in fruit firmness.

[IMAGE OMITTED. SEE PDF]

Development and validation of a KASP assay for prediction of fruit firmness

Significant SNPs associated with the QTL on 6A were tested for KASP designing and resulted in a subgenome-specific assay for Axiom SNP AX-184242253. We tested the predictive capacity of KASP-184242253 on a diverse collection of 138 strawberry accessions evaluated for fruit firmness during an additional fourth season (2022–2023; Table S8). TT, CT, and CC genotypic classes showed significant differences in fruit firmness (Figure 8). Strawberry accessions with the homozygous alternative allele (TT) produced fruits 48.13% softer than those with the reference CC allele. Heterozygous genotypes displayed an intermediate phenotype. It is noteworthy that only three accessions out of the 59 from subpopulation 1 had the TT genotype, clearly indicating a breeding selection on this locus. Similarly, subpopulations 3 + 4 + 5, encompassing 56 older accessions, only included two with the CC genotype. Comparison of genotype means within subpopulation 1 and within the subgroup including subpopulations 3 + 4 + 5 resulted in 34.06% and 42.19% difference in fruit firmness between TT and CC genotypes, respectively. A GWAS for fruit firmness using the 138 accessions and the 2022–2023 phenotypic data resulted once more in the detection of the same QTL on 6A (Figure S15). The phenotypic variance explained by this QTL was 44% and the allelic effect was −0.426 N. Together, these results indicate the great utility of this KASP assay for marker-assisted selection (MAS).

[IMAGE OMITTED. SEE PDF]

DISCUSSION

Analysis of genetic diversity and population structure showed a diversification into six subpopulations, with considerable admixture and minor structuring

Population structure and PCA indicated a diversification based on the breeding history and geographic adaptation, as in other reported analyses (Gil-Ariza et al., 2009; Zurn et al., 2022). Six subpopulations were identified, as in another recent study involving 1300 octoploid strawberry accessions of cultivated and wild species (Hardigan, Lorant et al., 2021). Older European varieties and the more recent Californian-Mediterranean type varieties (subpopulations 1 and 3 in Figure 1C) were clearly separated, while three intermediate subpopulations contained much more genetic admixture (subpopulations 4, 5, and 2 in Figure 1C). Both PCA and NJ tree separated the recent hybrids with F. chiloensis (Figure 1A,C). Despite these results, the low F_ST value suggests that the different subpopulations still share considerable genetic diversity, and thus the collection is not very structured, which reduces the probability of false positives in GWAS. The value obtained for our population was similar to that reported in other studies using other octoploid populations (Hardigan, Lorant et al., 2021; Zurn et al., 2022).

Relationship between traits and PCA highlights the main factors affecting phenotypic variation in the collection

In the three seasons, a large phenotypic variation was observed among individuals for all the evaluated traits, which was generally consistent among the 3 years (Table 1). According to this, positive correlations between seasons were in general medium-high. Correspondingly, phenotypic variation of many traits was also affected by different environmental factors between seasons. Correlations and PCA of phenotype data highlighted interesting relationships between traits (Figure 2; Table S3). First, the most vigorous genotypes produced the larger fruits, a higher number of fruits (therefore, greater yields), and also a higher number of runners. Similar conclusions were recently reported for another strawberry germplasm collection (Hummer et al., 2023). Our study also highlights a positive relationship between those traits and firmer fruits. In fact, fruit weight and firmness strongly impacted the phenotypic variation observed among accessions and across the strawberry breeding history (Figure 2B,C; Figure S1). Genomic signatures of domestication have been recently reported for these two breeding targets (Fan & Whitaker, 2023; Hardigan, Lorant et al., 2021). In agreement, the lead SNPs for fruit weight QTLs on chromosomes 3C and 3D (and not for those on 2B or 4C), and those for fruit firmness QTL on 6A (and not those on 3A) were within or less than 0.25 Mb from selective sweeps described in both Hardigan, Lorant et al. (2021) and Fan and Whitaker (2023).

Earlier runnering dates were correlated with higher production of runners. Also, more vigorous plants seemed to runner earlier. Negative correlations were also observed between runnering time and fruit number, indicating competition between vegetative and reproductive growth for available meristems (Gaston et al., 2013; Muñoz-Avila et al., 2022). Positive correlations between plant diameter and/or height with the number of runners have been previously reported (Antanaviciute et al., 2017; Hummer et al., 2023).

Fruit color traits were all significantly correlated and contributed greatly to the large diversity observed in the population (Figure 2; Table S3). As in other works (Ali & Serçe, 2022; Lerceteau-Köhler et al., 2012; Zorrilla-Fontanesi et al., 2011), we observed a strong positive correlation between fruit external lightness and yellow external hue. Both the external color and the internal color of the pulp (including the core) were positively correlated, which has been reported in other works, and is to be expected since all three depend on an accumulation of anthocyanins (Castillejo et al., 2020).

Strawberry accessions with higher SSC (FR_Brix) had fruits with higher vitamin C content (Figure 2; Table S3) in agreement with other reports (Ali & Serçe, 2022; Zorrilla-Fontanesi et al., 2011). Strikingly, these accessions belong mostly to the old European and American varieties (subpopulation 3), suggesting that breeding had a negative impact on these two important traits. In agreement with this conclusion, SSC and AsA were negatively correlated with fruit number, weight, and firmness, the main breeding targets throughout history. Negative correlations between these traits have also been reported in other strawberry populations (Cockerton et al., 2021; Hummer et al., 2023; Whitaker et al., 2012; Zorrilla-Fontanesi et al., 2011). As previously reported, positive correlations were observed also between SSC and fruit acidity (Vallarino et al., 2019; Wada et al., 2020). This positive correlation could in part be caused by the fact that, besides sugars, pectins, organic acids, and amino acids also contribute to Brix values.

Genotype-phenotype associations varied depending on the season and the GWAS model

Our GWAS for agronomic and fruit quality traits has detected 121 significant SNPs associated with 19 of the 26 evaluated traits. Many of them represent novel and distinct QTL, while others overlap previously identified regions. Previously identified regions, as for fruit color traits on 1B caused by natural variation on FaMYB10 (Castillejo et al., 2020), validate the experimental methods used, while novel associations, particularly those detected in more than one season and/or with larger phenotypic effects, represent the starting point for the identification of the underlying genes. Also, significant markers with large effects here detected for many key breeding traits, such as runnering date and runner number, fruit weight, or firmness, represent useful tools for MAS.

Association mapping is characterized by high resolution and power to detect monogenic traits, such as many resistance genes (Bush & Moore, 2012). However, for polygenic traits where many loci explain small amounts of phenotypic variation, the statistical signals are much weaker and more difficult to detect than for genetically simpler traits (Atwell et al., 2010). This might be the reason why no significant associations were detected for complex polygenic traits such as plant height, flower number, flowering time, or AsA content. Studying these complex traits would require a higher number of accessions (Balding, 2006; Bush & Moore, 2012), which may have limited the statistical power with our experimental population size. Twelve of the traits were evaluated using scores, including FL_GCP, which was evaluated with a binomial scale. Although the use of scales might reduce the power of the association tests (Poland & Nelson, 2011), only in four of them, including FL_GCP, no significant associations were detected. However, petal greening was the only trait without significant correlations between seasons, indicating that it might be mainly controlled by environmental factors.

Despite reasonable correlations between seasons and a similar range of variation across years (Table 1), most associations were detected in individual seasons, suggesting substantial environmental effects. GWAS results are not always consistent between years, as those reported by McClure et al. (2018) for fruit quality traits in apple (Malus domestica Borkh.), in which fruit firmness and color displayed high correlations between years, but different associations were detected for each year. In our study, most of the associations were detected in the third season (Table 3; Table S4), which may have been influenced by the larger sample size. The use of four models allowed the evaluation of the best fit to our population, trying to avoid possible interactions that could produce false positives, or the lack of significant associations by false negatives. The BLINK model resulted in the highest number of significant SNPs and, in general, produced QQ plots with the best connection between observed and expected p-values, indicating an effective control of population structure and family relatedness (Table 3; Table S4).

Differential expression of FaPG1 controls natural variation in fruit firmness in strawberry

Improvement of fruit firmness has been a key breeding target since the origin of the species, about 300 years ago, due to its effects both on the organoleptic quality (texture) and postharvest and shelf life of the fruit. For this trait, we identified seven QTLs, although two of them, on chromosomes 3A and 6A, had larger phenotypic effects (Figure 6; Table S4). The QTL on chromosome 6A was detected in the third season, the one with the larger accession sample, using the four models and increased fruit firmness by about 0.27 N. Significant associations in exactly the same region have been detected in other populations (Cockerton et al., 2021; Hardigan, Lorant et al., 2021). We further detected this QTL using a larger sample size and the FaRR1 genome in an additional fourth season (Figure S15). The phenotypic variance explained (44%) and its detection in other populations and environments, indicate the usefulness of associated SNPs for MAS.

The 891-kb long haploblock enclosing the significant markers contains the FaPG1 gene. Downregulation of FaPG1 expression in the cultivar ‘Chandler’ by antisense reduced strawberry fruit softening at ripe and postharvest stages (Garcia-Gago et al., 2009; Quesada et al., 2009), highlighting its key role on this trait. Further studies in FaPG1-antisense fruits revealed a 42% decrease in pectin solubilization and reduced depolymerization of tightly bound pectins (Posé et al., 2013). More recently, FaPG1 knockout mutants generated by CRISPR/Cas9 editing resulted in fruits with increased firmness, reduced softening rate during postharvest, and less susceptibility to Botritis cinerea damage (López-Casado et al., 2023). Genotyping using the KASP-184242253 assay indicated that while ‘Camarosa’ was homozygous for the C allele, ‘Chandler’ was heterozygous and, in agreement, displayed softer fruits (Table S8). We have shown a 13-fold change in FaPG1 expression between strawberry accessions contrasting in fruit firmness (Figure 7). We have also shown significant differences in fruit firmness across the six subpopulations of the germplasm collection (Figure S1), which agrees with the identification of the same QTL on 6A overlapping with early-phase domestication sweeps (Hardigan, Lorant et al., 2021). Finally, the KASP assay here developed was able to predict a 48.13% increase in fruit firmness. This effect was still larger than 34% when analyzing modern or older accessions independently. As in previous reports, our study highlighted the higher SSC and lower firmness in older cultivars. The KASP-184242253 assay represents an extremely useful tool for MAS, particularly in breeding programs aiming to increase fruit sweetness without detriments to firmness. Together with previous knockdown experiments (Garcia-Gago et al., 2009; Quesada et al., 2009), our results indicate that FaPG1 is the main contributor to natural variation in strawberry fruit firmness.

CONCLUSION

Breeding superior fruit cultivars involves the simultaneous selection of many agronomic and fruit quality traits. In this study, the genetic control of 26 traits was studied by GWAS using a diverse collection of 124 strawberry accessions. Our results highlighted negative relationships between fruit weight and firmness, with SSC, acidity, and AsA over strawberry breeding history. Genotype–phenotype association studies identified a total of 95 QTLs for key breeding traits including runnering, fruit weight, color, SSC, acidity, and firmness. Transcriptional regulation of FaPG1 accounted for 47% natural variation on fruit firmness and underlies the major QTL for fruit firmness on chromosome 6A. While down-regulation of FaPG1 by biotechnological means increases fruit firmness and extends postharvest life (López-Casado et al., 2023), the KASP marker linked to FaPG1 here developed represents an extremely useful alternative tool for MAS. Overall, the findings provide valuable insights into the genetic control of agronomic and fruit quality traits and provide associated markers for molecular breeding of superior strawberry cultivars that are more productive and with improved fruit quality.

AUTHOR CONTRIBUTIONS

Pilar Muñoz: Data curation; formal analysis; investigation; visualization; writing—original draft. Francisco Javier Roldán-Guerra: Data curation; formal analysis; investigation; visualization. Sujeet Verma: Software. Mario Ruiz-Velázquez: Investigation. Rocio Torreblanca: Investigation. Nicolás Oiza: Investigation. Cristina Castillejo: Formal analysis; methodology; supervision. José F. Sánchez-Sevilla: Conceptualization; data curation; formal analysis; funding acquisition; software; supervision. Iraida Amaya: Conceptualization; formal analysis; funding acquisition; supervision; visualization; writing—review and editing.

ACKNOWLEDGMENTS

This work was supported by Ministerio de Ciencia e Innovación and Agencia Estatal de Investigación (PID2019-111496RR-I00 and PID2022-138290OR-I00/MCIN/AEI/10.13039/501100011033/FEDER), Junta de Andalucía and FEDER (P18-RT-4856) and by the European Union's Horizon2020 research and innovation program (BreedingValue project, grant agreement 101000747). The Fragaria collection at IFAPA is financed by IFAPA Project PR.CRF.CRF202200.002 with funds from the European Agricultural Fund for Rural Development. Pilar Muñoz acknowledged a PhD FPI-INIA contract co-financed by FSE. We are grateful to Francisco J. Durán for his excellent care of strawberry plants.

CONFLICT OF INTEREST STATEMENT

The authors declare no conflicts of interest.

DATA AVAILABILITY STATEMENT

All data generated or analyzed during this study are included in this published article (and its supplementary information files).

References

Ali, M. N., & Serçe, S. (2022). Vitamin C and fruit quality consensus in breeding elite European strawberry under multiple interactions of environment. Molecular Biology Reports, 49, 11573–11586. [DOI: https://dx.doi.org/10.1007/s11033-022-07849-5]