A pseudo full-sibling F₁ population of 165 individuals obtained from a cross between the variety “Capitola” (“CA75.121-101” × “Parker,” University of California, Davis, CA, United States) and the advanced line “CF1116” {[“Pajaro” × (“Earlyglow” × “Chandler”)], reference from the Ciref, France} was developed. The “Capitola” and “CF1116” parents display contrasting fruit colour together with differences in fruit shape and weight, firmness, sweetness, and acidity (Lerceteau-Köhler et al., 2012). For each of the two consecutive study years (2010 and 2011), cold-stored strawberry plants planted in 2009 and 2010 were grown in soil-free pine bark substrate under plastic tunnel with daily ferti-irrigation and control of biotic stresses. The mapping population included a total of 165 individuals over the two study years. Within this progeny, 72 and 131 individuals, including the parents, were respectively phenotyped in 2010 and in 2011. Fruits were harvested at the red ripe stage, when red colouration of the fruit is homogeneous, and processed as previously indicated (Labadie et al., 2020) to produce frozen powder samples that were further stored at −80°C until use for chemical analyses.

In 2011, for both harvests, red ripe fruits (4--5 fruits per harvest) were also photographed side-by-side with the strawberry colour chart from Ctifl.¹ Fruit colour was then scored by two independent persons (two replicates) on a scale from 0 (very pale red-orange) to 6 (very dark red) using the photographs. A mean score value was then obtained for each genotype (Supplementary Figure 1).

Analysis of polyphenolic metabolites by LC–ESI-MS was done as previously described (Ring et al., 2013). A total of 16 traits encompassing 13 individual phenolic metabolites and their sum by chemical class (three traits) were measured for the 2 years as described in Labadie et al. (2020). Analyses of pooled frozen powder samples were carried out in 2010 and 2011 on six replicates for the parents and on three replicates for individuals from the progeny.

Extraction and measurement of total anthocyanin content (ANTHc) by colourimetric assay were as described in Labadie et al. (2020). Results are expressed as mg PgGs equivalents/100 g fresh weight. For each genotype, four technical repeats from the pooled two-harvest-fruit-powder were performed.

Frequency distribution of each trait was represented using ggplot2 (Wickham, 2016) (v3.2.1) r-package.

DNA from the parental lines “Capitola” and “CF1116” and from 165 individuals from the mapping population was extracted using DNeasy Plant Mini kit (Qiagen, Hilden, Germany) and hybridised with the Affymetrix^® 90 K Axiom^® SNP array (Affymetrix, CA, United States) at CeGen USC (Santiago de Compostela, Spain). IStraw90^® SNP array includes 138,000 SNP markers corresponding to 90,000 localisations on wild diploid strawberry F. vesca genome (Bassil et al., 2015; Edger et al., 2020). Analysis was performed using Genotyping console™ and SNPpolisher© (Affymetrix, CA, United States) following manufacturers recommendations.

We developed a molecular marker, named BxANR_5UTR, in the 5′UTR of the ANR gene. To this end, using the F. vesca reference genome and WGS sequencing data from both parents, “Capitola” and “CF1116,” we designed primers for BxANR_5UTR (forward: GAGAGTTGGTGGTGCTTTCA; reverse ATGGT GTGGGTGTGTCTCAG, 5′-3′) to amplify a 182 bp fragment that spans an 18 bp gap (TTCTTCCTCTTCTTCTTC) in the ANR 5′UTR. As described in Perrotte et al. (2016), we used specific extended primers of the M13 sequence which allow the hybridisation of a fluorolabelled primer (6-FAM or VIC). The PCR products of all individuals of the segregating population were diluted 1:200 and separated using capillary electrophoresis [CE; Applied Biosystems (ABI 3730), ABI Genescan™ 500 LIZ size standard]. Peak identification and fragment sizing of each electropherogram were done with GeneMapper v4.0 software. Data were converted to a qualitative binary code for scoring the presence or absence of BxANR_5UTR homoeo-alleles.

Single dose markers (SD) from the Affymetrix array (Bassil et al., 2015) that were in backcross configuration and segregated 1:1 (Rousseau-Gueutin et al., 2008) were used in combination with previously mapped SSR, SSCP, and AFLP markers (Gaston et al., 2013; Labadie et al., 2020) for map construction using JoinMap 5.1 software (Van Ooijen, 2011). Grouping was performed using independence log of the odds (LOD) and the default settings in JoinMap^®. LGs were chosen from an LOD higher than 10 for all of them. Map construction was performed using the maximum likelihood (ML) mapping algorithm and the following parameters: chain length 5.000, initial acceptance probability 0.250, cooling control parameter 0.001, stop after 30.000 chains without improvement, length of burn-in chain 10.000, number of Monte Carlo EM cycles 4, chain length per Monte Carlo EM cycle 2.000 and sampling period for recombination frequency matrix samples: 5.

For QTL analysis, the female and male linkage parental maps based on the 165 individuals were used separately. Phenotypic data of the 72 and 131 individuals for 2010 and 2011, respectively, were represented by the mean value of the three replicates. QTL detection was performed by simple interval mapping (SIM) using R/QTL (Broman et al., 2003). Permutation analysis (1,000 permutations) was performed to calculate the critical LOD score. QTL with LOD values higher than the LOD threshold at P ≤ 0.05 were considered significant. When one QTL was found significant, we used composite interval mapping (CIM) with one co-variable at the position of the significant QTL and reiterated the analysis until no new significant QTLs were detected. Bayesian credible interval was calculated using the function “bayesint” at probability of 0.95. The proportion of phenotypic variance explained by a single QTL was calculated as the square of the partial correlation coefficient (R²). Mapping results are displayed using MapChart (Voorrips, 2002).

A custom-made oligonucleotide-based (60-mer length) platform (Roche NimbleGen) designed from non-redundant F. vesca strawberry sequences (Fv_v1.0, Shulaev et al., 2011) representing a total of 18,152 unigenes (Ring et al., 2013) was used for the analysis of differentially expressed genes between 21 F1 genotypes from the “Capitola” × “CF1116” population displaying contrasted phenotypic values for flavonoids and colour-related traits. RNA extraction from red ripe fruits and microarray processing and data analysis were as previously described (Ring et al., 2013). Student’s t-test was used with a confidence of P < 0.05 to detect statistically significant differences.

Whole genome sequencing of the two parents “Capitola” and “CF1116” was performed using paired-end Illumina sequencing with a ∼ 50X coverage of the F. × ananassa genome. Illumina paired-end shotgun-indexed libraries were prepared and sequenced using an Illumina HiSeq 3000 at the Institut National de la Recherche Agronomique GeT-PlaGe facility (Toulouse, France), operating in a 150-bp paired-end run mode. Raw fastq files were mapped to the F. vesca FvH4 reference genome sequence (F. vesca Genome v4.0.a1) (Edger et al., 2018) using BWA-MEM algorithm (Li, 2013) for the alignment of paired-end (150 bp) Illumina reads. Polymorphisms between “Capitola” and “CF1116” were identified using Integrative Genomics Viewer (IGV) (Robinson et al., 2011). All identified polymorphisms were tested for the 1:7 (mutant allele:WT alleles) segregation ratio for goodness-of-fit to theoretical ratio (Chi-squared test) when considering the hypothesis that, out of the eight homoeo-alleles expected in the octoploid F. × ananassa, one single homoeo-allele (mutant allele) specific to “Capitola” or to “CF1116” controls the trait.

In addition to the thirteen fruit flavonoid compounds previously analysed by LC-ESI-MS (Labadie et al., 2020; Figure 1A), we mapped the colour of the fruit (COLOUR trait) scored visually according to the strawberry colour chart (Ctifl² (Supplementary Figure 1). As described in Labadie et al. (2020), the flavonoids detected in the fruit were flavonols (four compounds), flavan-3-ols (four compounds), and anthocyanins [five compounds: PgGs, pelargonidin-3-glucoside-malonate (PgGsM), pelargonidin-3-rutinoside (PgRs), cyanidin-3-glucoside (CyGs), and (epi)afzelechin-pelargonidin-glucoside (EpPgGs)]. The quantitative value and range in parents and progeny, the broad sense heritability and the transgression value in the “Capitola” × “CF1116” population were reported in Labadie et al. (2020) for all the traits except COLOUR that was measured only in 2011 (Supplementary Table 1 and Supplementary Figure 1). Very high broad sense heritability values ranging from 0.81 (PgGs) to 0.93 (PgGsM) were observed for the major anthocyanin compounds (Labadie et al., 2020). Distributions of flavonoids and colour-related traits are shown in Figure 1B for 2010 and in Supplementary Figure 2 for 2011. Continuous variations of the phenotypic values were observed in the progeny for all the traits assessed in 2010 and in 2011.

In this study, we focused on colour-related mQTLs and QTLs. The PgGs, which contributes to total ANTHc for as much as ∼90%, displayed a fivefold variation among the most contrasted individuals. This is in agreement with the distribution of most flavonoid metabolites and colour-related traits in the progeny, which showed a 4–10-fold variation. As could be expected, the ANTHc values measured by colourimetry varied in the same range as PgGs. In contrast, considerable variations were observed for the minor anthocyanin forms CyGs and PgRs, which displayed 23-fold and 114-fold variations in the progeny, respectively; and for the COLOUR score, which ranged from 0.5 to 6. PgGs -malonate and CyGs were not detected in a large number of individuals, resulting in a very skewed distribution in the progeny (Figure 1B and Supplementary Figure 2).

For linkage map construction, we added to the previous linkage maps (Rousseau-Gueutin et al., 2008) 9,455 SNP markers from the Axiom^® IStraw90^® SNP array [only single dose markers (SD) with 4,974 female and 4,481 male markers] (Bassil et al., 2015). The construction of the linkage maps was done with a total of 5,216 and 4,879 markers for the female and male linkage maps, respectively. The final number of markers covered the expected 28 LGs for the female map and 31 LGs for the male map (3 additional small LGs were included) (Supplementary Tables 2–4). The lengths of the female and male linkage maps were 4,135.1 and 3,929.2 cM, respectively, with an average distance between markers of 0.8 cM.

Linkage groups were assigned to one of the seven homoeologous groups (HGs). However, contrary to our previous reports (Rousseau-Gueutin et al., 2008; Labadie et al., 2020), a homoeologous LG (e.g., LG1A, B, C, or D) has been named according to the recent nomenclature of Hardigan et al. (2020, 2021a) where letters refer to species-derived subgenomes: A, F. vesca; B, F. iinumae; C, F. nipponica, and D, F. viridis (Supplementary Table 2).

QTL were detected for all the quantitative traits analysed using CIM for “Capitola” or “CF1116.” Information on markers for the male and female maps are provided in Supplementary Tables 3, 4. The list of significant QTLs detected for each trait including associated markers, position on the male and female linkage maps, LOD score and effect is provided in Supplementary Table 5. The QTL values at the 5 and 10% thresholds used to select the significant male and female QTLs for each trait and year are given in Supplementary Table 6. Distributions per HGs and LGs of significant colour-related mQTLs [PgGs, PgGsM, PgRs, CyGs, EpPgGs, and Ant (summation of anthocyanins compounds)] and QTLs (ANTHc and COLOUR) are synthesised in Table 1 for the male and female. The graphical representation of the location on the male and female linkage maps for the significant QTLs detected for the 2 years of study (2010 and 2011) is presented in Supplementary Figure 3.

Among the 65 QTLs detected for all the flavonoid and colour-related traits analysed (Supplementary Table 5), a total of 37 (57%) significant colour-related QTLs (29 anthocyanins mQTLs, 4 ANTHc, and 4 COLOUR QTLs) was detected (Tables 1, 2). More than half of the colour-related QTLs (54%) were located on the F. vesca subgenome. Most LGs harboured only one region controlling variations in colour-related traits, the notable exception being the LG6A carrying two male and three female regions with colour-related QTLs (Table 1). Quantitative variations of few anthocyanin compounds were controlled by mQTLs located on different homoeologous LGs within the same HG (Table 2). The only examples are those of PgGs QTLs located on the female LGs F3A and F3D in 2011 and of PgRs QTLs located on the male LGs M5A and M5C and M7B and M7C in 2011. Moreover, half of the colour-related QTLs detected on the male linkage map were located on only three LGs (M1B, M3A, and M6A) (11 QTLs out of 22). On the female linkage map, the overwhelming majority (12 QTLs out of 15) were located on only three LGs (F2A, F3A, and F6A). Remarkably, the above-mentioned LGs correspond to the F. vesca subgenome. Analysis of the contribution of the different subgenomes further highlighted the predominant role of F. vesca in the determination of fruit colour in F. × ananassa. The F. vesca subgenome is responsible for major variations in PgGs content (F2A QTL: R² = 15; M3A QTL: R² = 20.6; F3A QTL: R² = 8.8; F6A QTL: R² = 19.7) and contributes to 42% of the total number of colour-related QTLs detected (Table 1).

The F. iinumae (17%), F. nipponica (17%), and F. viridis (21%) subgenomes each contribute less than 25% of all 37 colour-related QTLs (Table 1). Subgenomes derived from F. iinumae (M1B QTL; R² = 11.65) and from F. viridis (M3D and M4D QTLs; R² > 7.5) contribute to variations of PgGs content in the progeny (Table 2 and Supplementary Table 5). F. iinumae, F. nipponica, and F. viridis are also involved in the control of PgRs and F. viridis in that of COLOUR (F4D; R² = 9). We further investigated if the MYB10-2 homoeo-allele located on the LG1 F. iinumae–derived subgenome (LG1-2, i.e., LG1B) (Castillejo et al., 2020) could underline the M1B Ant, PgGs, PgRs, AfPgGs, ANTHc, and COLOUR QTLs (Supplementary Table 5). In the male linkage map, the closest marker to the MYB10-2 homoeo-allele at position 15,517,937 bp on the subgenome Fvb1-2 of the Camarosa reference genome is AX-89846847 at 118.78 cM on M1B. This marker is in the Bayesian credible interval of the flavan-3-ols AfCat_2011 mQTL (Supplementary Figure 3) but not in the other M1B QTLs intervals detected in our study.

The two colour-related QTLs showing the highest percentage of variance (more than 20%), which are located on M3A and F6A (Ant_2010, PgGs_2010, and COLOUR_2011) (Table 2), display opposite effects on ANTHc (Supplementary Table 5). We investigated if they had epistatic relationships through the analysis of the allelic interaction by ANOVA (Figure 2). Variance analysis confirmed that the male M3A QTLs had indeed a positive effect on trait values and the female F6A QTLs a negative effect. However, because trait values were similar in individuals harbouring both M3A and F6A QTLs, we can conclude that there is only additive effect and likely no significant interaction between them.

We also identified overlapping Bayesian credible intervals for various QTLs, which would indicate the presence of pleiotropic or closely linked QTLs (Table 2 and Supplementary Figure 3). In particular, overlapping QTLS were found for total anthocyanins (Ant), PgGs and COLOUR on M3A and F3A (Figure 3A and Supplementary Table 5). We could detect mQTLs for Ant and PgGs on the male map (M3A) in 2010 and on the female map (F3A) in 2011, and a COLOUR QTL on M3A in 2011 (Figure 3A). CIM analysis further showed that the fruit COLOUR QTL and Ant and PgGs mQTls, which displayed high LOD scores values >3.0, were co-located in a narrow chromosomal interval (Figures 3B–D). Furthermore, the male allele from “CF1116” (on M3A) had a positive effect on the levels of anthocyanins and COLOUR (Figures 3B,C) while the female allele from “Capitola” (on F3A) had the opposite effect (Figure 3D). This points to the likely presence in this LG3A region of two different allelic variants, one from the male parent and one from the female parent, which affect either positively (the male homoeo-allelic variant) or negatively (the female homoeo-allelic variant) the colour-related traits. For Ant and PgGs, the percentages of phenotypic variance explained by the M3A mQTLs (R² = ∼20) were more than twice as high as that of the F3A mQTLs (R² = ∼8.5) (Supplementary Table 5).

To detect candidate genes underlying colour QTLs and their homoeo-allelic variants, we took advantage of the strong synteny between F. vesca genome (Edger et al., 2018) and the subgenome of F. vesca Fvb3-4, which corresponds to our LG LG3A (Hardigan et al., 2021a). For both QTLs, we identified on the diploid FvH4 F. vesca genome an interval framed by SNP markers with physical positions overlapping the Bayesian credible intervals of the M3A and F3A QTLs. On M3A, this region is flanked by AX-89904962 and AX-89785774 Affymetrix markers and spans an interval on chromosome 3 (Fvb3) from 1,213,489 to 2,673,762 b, while on F3A this region is larger (826,085–2,539,615 b) and almost overlapped with the M3A interval. Based on the latest annotation of the F. vesca genome (Li et al., 2019), we identified a total of 392 genes in the M3A/F3A interval (826,085–2,673,762 b) (Supplementary Table 7) and searched them for candidate genes possibly involved in the regulation and/or synthesis of flavonoids in strawberry fruit. Among them, the most likely candidate gene is annotated as NAD(P)-binding Rossmann-fold superfamily protein (FvH4_3g02980). It encodes an ANR enzyme, which catalyses the conversion of pelargonidin to epiafzelechin and of cyanidin to epicatechin (Figure 1A) and has previously been demonstrated to control anthocyanin accumulation in F. vesca fruit (Fischer et al., 2014). Two additional candidate genes encoding MYB-related transcription factors annotated as MYB58 (FvH4_3g03680) and MYB102-like ODORANT (FvH4_3g03780) are found in the M3A/F3A interval. To further investigate the candidate genes, we first mined transcriptome data obtained by microarray analysis of 21 individuals from the segregating population which displayed contrasted flavonoid-related phenotypes. To this end, we first constituted two groups of seven F1 individuals each displaying either high or low PgGs content at the red ripe stage, with a PgGs content ratio between the two pools of 1.7 in 2010 and 1.4 in 2011 (Figure 4 and Supplementary Table 8). We next analysed these individuals for differentially-expressed genes (DEGs) using a custom-made oligonucleotide-based (60-mer length) platform representing a total of 18,152 strawberry unigenes (Ring et al., 2013). To identify significant DEGs present in the region of interest, we applied a Student’s t-test on the two phenotypic groups for all genes located in the Bayesian credible interval. Out of the 392 genes found in the M3A/F3A interval that harbours the QTL and mQTLs of interest, 304 genes including ANR (FvH4_3g02980), MYB102-like ODORANT (FvH4_3g03780) and MYB58 (FvH4_3g03680) were present on the microarray (Supplementary Table 9). Among them, 50 DEGs (P-value < 0.05) were found, which included ANR and MYB102-like ODORANT but not MYB58. The expression of the ANR and MYB102-like ODORANT genes were respectively 1.6 and 1.9-fold higher in the individuals with low PgGs content than in the individuals with high PgGs content (Figure 4 and Supplementary Table 9). Analysis of the remaining DEGs did not highlight any additional candidate gene for the control of anthocyanin accumulation in the fruit (Supplementary Table 9).

We then performed the whole genome sequencing (Illumina Hiseq 3000) of the “Capitola” and “CF1116” parents of the segregating population. We obtained 137 and 142 million 150 pb paired-end reads, representing 50.7X and 52.7X coverage of the octoploid genome (Edger et al., 2019), respectively. As summarised in Figure 5A, the LG3A colour QTL marker has a positive effect on fruit colour in the male (“CF1116”) and a negative effect in the female (“Capitola”). Thanks to the alignment of paired reads “Capitola” and “CF1116” to the FvH4 F. vesca genome (Edger et al., 2018), which is syntenic to the subgenome Fvb3-4 of the Camarosa octoploid reference genome (Hardigan et al., 2021a), we had access to sequence polymorphisms of ANR and MYB102-like ODORANT genes in the four sub-genomes of both parents (Figure 5B). Because ANR and MYB102-like ODORANT are up-regulated in progeny individuals with low PgGs content (Figure 4), we first analysed sequence polymorphisms in regions that may affect gene regulation and that were different between both parents. In the 1 kb region upstream of the transcription start site of both genes, we found several SNPs and insertions by comparison of the “Capitola” and “CF1116” sequences. In the 5′untranslated region (UTR) of ANR, we found two deletions, one of 3 bp in Capitola and one of 18 pb in CF1116. In the 5′UTR of MYB102-like ODORANT, we found only one deletion of 8 bp in “Capitola.” We next analysed the protein coding regions for polymorphisms that are different between both parents and detected two SNPs in ANR and three SNPs in MYB102-like ODORANT.

To further reduce the number of candidate polymorphisms, we subsequently considered only those polymorphisms fitting the 1:7 ratio expected for single male or female homoeo-alleles (Chi-squared test; P > 0.05) located on a single LG (Figure 5C). Interestingly, the 18 bp deletion located from −109 to −127 bp upstream of the start codon in the 5′UTR of the ANR gene of “CF1116” (the male parent) (Figure 5C) was present in 13% of the reads. Analysis of the corresponding CTTCTTCCTCTTCTTCTT sequence with PlantRegMap (Jin et al., 2017) predicted that the deletion was in the 21 nucleotides binding site of a MADS Agamous-like transcription factor (Arabidopsis AT2G45660). Blast analysis of F. vesca v1.0 ab hybrid reference genome sequence at GDR (Jung et al., 2019) using AT2G45660 protein sequence as a query and visualisation of gene expression of top hits using F. vesca eFP browser (Hawkins et al., 2017) allowed the identification of four strawberry homologs expressed in fruit cortex along development (gene24852, gene04229, gene26119, and gene06301). Other small insertion/deletions (INDELS) or SNPs found in promoter or 5′UTR of ANR and MYB102-like ODORANT genes did not match known motifs. Two SNPs found in the protein coding region of ANR are synonymous and do not affect the function of the protein. One SNP found in the protein coding region of MYB102-like ODORANT is synonymous while two SNPs found in exon 3 of the MYB102-like ODORANT of “Capitola” (the female parent) are non-synonymous mutations leading to H146Q and I206N amino acid substitutions (Figure 5C).

To develop genetic markers for improvement of strawberry fruit colour by MAS, we first analysed the allelic status of SNP marker AX-89826853 (position: 1,981,447 bp on FvH4, i.e., close to ANR) using the 14 selected progeny individuals with contrasted PgGs content. All individuals with high PgGs showed the presence of this SNP marker while individuals with low PgGs did not. To further check whether the 5′UTR deletion found in ANR could be used for developing a genetic marker tightly linked to the M3A colour QTLs, we next PCR-amplified a 182 bp region (size in the diploid reference genome) spanning the 18 bp deletion found in “CF1116” and analysed the resulting PCR product by CE (Figure 6A) in the segregating population derived from “Capitola” and “CF1116.” Beforehand the PCR product (hereafter named BxANR_5UTR) was cloned in a plasmid vector and sequenced to confirm that the expected 5′UTR sequence of ANR was amplified. We identified five CE peaks in “Capitola” and “CF1116” that correspond to five homoeo-alleles; three of them are common to both parents (h, k, and e), one is specific to “Capitola” (f) and one is specific to “CF1116” (g) (Figure 6B). Noteworthy, the estimated size of the peak for allele g, which was observed only in the male parent “CF1116,” corresponds approximately to the size of the deleted variant of BxANR_5UTR while the size of the allele f would correspond to the 3 bp deleted variant specific to “Capitola” (Figure 6A). The h and k homoeo-alleles found in both “Capitola” and “CF1116” have smaller peak sizes corresponding to PCR fragments carrying several deletions.

Within the five homoeo-alleles, three homoeo-alleles (e, f, and g) and two homoeo-alleles (h and k) are in coupling/repulsion configurations (ef × eg and hk × hk). Using JoinMap, we further mapped the homoeo-alleles to LG3A (Fvb3-4; alleles e, f, and g) and LG3B (Fvb3-2; alleles h and k) (Figure 6B). Surprisingly, BLAST search of the “Camarosa” reference genome sequence (Jung et al., 2019) found only one ANR gene, which is located on Fvb3-2, and failed to identify ANR on Fvb3-4, which is the subgenome harbouring the M3A colour QTL. Thanks to the position of Affymetrix markers close to the e, f, g, h, and k homoeo-alleles, we could nevertheless confirm that the positions of the e, f, and g homoeo-alleles are orthologous to those of the h and k homoeo-alleles, thus indicating that the ANR gene is indeed present in the Fvb3-4 subgenome of “Capitola” and “CF1116.” Possible explanations are that ANR has been deleted from Fvb3-4 in the “Camarosa” variety or that the sequence of the Fvb3-4 region is incomplete in the reference genome. We favour the last hypothesis because this gene is present in the chr_3A (corresponding to Fvb3_4) of the recently released genome sequence of Royal Royce (Hardigan et al., 2021b).

We next analysed the effect of the various combinations of the five homoeo-alleles on fruit colour intensity in the segregating population derived from “Capitola” and “CF1116” (137 progeny individuals and the two parents) (Figures 6C,D). Remarkably, the g homoeo-allele, which corresponds to the deleted variant of BxANR_5UTR that is co-localised on M3A (Fvb3-4) with the AX-89826853 marker linked to the PgGs, Anth and COLOUR QTLs (Figures 3A,B, 6B), had a significant and positive effect on fruit colour as expected (Figure 5A). When present in the genetic combination, the g homoeo-allele resulted in darker red fruit with an average increase in the COLOUR score from 3 to 3.9 (on a 0 to 6 scale; Supplementary Figure 1), which is considerable (Figures 6C,D). Because the analysis of BxANR_5UTR is straightforward and allows the identification of strawberry genotypes with dark red fruit, it can therefore be used in MAS as a genetic marker for breeding strawberry varieties with more intense red fruit colour.

In polyploid plant species, each trait is likely controlled by homoeologous gene series or homoeo-alleles. In the octoploid strawberry, homoeo-alleles are located at orthologous positions on one of the four subgenomes (Edger et al., 2019). Remarkably, our results showed that, within a given HG, mQTLs for flavonoids and anthocyanins were mostly located on a single LG (i.e., on a single subgenome). As an example, the A subgenome was likely responsible for the major PgGs mQTLs localised on M3A, F3A, and F6A LGs. Furthermore, the main LGs accounting for half of male colour-related QTLs (three LGs) and for the majority of female QTLs (three LGs) can be attributed to the A subgenome, which is mainly derived from F. vesca (Hardigan et al., 2021a), thus highlighting the predominant role of F. vesca in fruit flavonoid metabolism and quality. Our findings are thus consistent with transcriptome analyses showing that F. vesca homoeologs are responsible for almost 89% of the anthocyanin biosynthesis in octoploid strawberry (Edger et al., 2019).

The contribution of the various diploid genomes to octoploid genome is complex (Edger et al., 2019; Liston et al., 2020, Hardigan et al., 2021a). We followed Hardigan et al. (2021a) for the assignment of a subgenome to a given LG. In the population studied, subgenomes other than F. vesca may play key roles in flavonoid metabolism, including, e.g., subgenomes derived from F. iinumae and F. viridis for PgGs; subgenomes derived from F. iinumae, F. viridis, and F. nipponica for PgRs; subgenomes derived from F. iinumae- and F. viridis for COLOUR. The crucial role in the regulation of strawberry fruit colour of the MYB10-2 homoeo-allele localised on LG1B (F. iinumae–derived subgenome) has been demonstrated recently (Castillejo et al., 2020). However, in our population, this homoeo-allele was not included in the Bayesian credible interval of the colour-related QTLs detected on M1B (Supplementary Figure 3) and is likely not responsible for the colour variations observed.

In F. × ananassa, pinpointing with high resolution the position of flavonoid mQTLs and other colour-associated QTLs opens the possibility to identify the homoeo-alleles responsible for major variations in colour-related traits. Major QTLs detected on LG3A for the total anthocyanins, the PgGs and the visually scored fruit colour are all co-located at the beginning of LG3A (Figure 3A) and overlap with the colour-related QTLs previously detected over 3 years of study for the physical parameters L and b (colour space values) (Lerceteau-Köhler et al., 2012), indicating that these colour QTLs are robust. We could further map with high resolution the genetic architecture of the colour trait and show that: (i) PgGs mQTLs displaying high LOD scores values >3.0 are co-localised on the male (M3A) and female (F3A) maps in a narrow chromosomal interval encompassing 392 genes, (ii) the male homoeo-allele from “CF1116” has a positive effect on PgGs content (Figures 3B,C) while the female homoeo-allele from “Capitola” has the opposite effect (Figure 3D), and (iii) homoeo-allelic variants of two candidate genes encoding ANR and MYB102-like ODORANT, both of which are differentially expressed in pools of individuals with contrasted PgGs contents, are respectively found in the male (“CF1116”) and female (“Capitola”) parents.

The possible candidate gene underlying the female F3A colour-related QTL (∼8.5% of explained variance) is a MYB transcription factor (TF). MYB TFs play predominant roles in the control of phenylpropanoid pathway and anthocyanin biosynthesis in strawberry (Salvatierra et al., 2013; Schaart et al., 2013; Medina-Puche et al., 2014; Wang et al., 2020). However, the most likely MYB candidate gene underlying the female F3A PgGs mQTL is not homologous to a well characterised strawberry MYB such as MYB10 (Medina-Puche et al., 2014; Hawkins et al., 2016; Castillejo et al., 2020; Wang et al., 2020; Zhang et al., 2020). It encodes a MYB102-like ODORANT protein homologous to the petunia R2R3-MYB ODORANT1 (ODO1) which regulates floral-scent related genes in petunia (Spitzer-Rimon et al., 2010). Ectopic expression of ODO1 in tomato activates phenylpropanoid metabolism without affecting volatiles (Dal Cin et al., 2011). In addition to the two SNPs leading to synonymous amino acid changes in ODO1 coding region, the MYB102-like ODORANT gene from “Capitola” LG3A carries several SNPs and INDELS in the 5′UTR and promoter regions. This opens the possibility that the differential expression of this gene, which was observed in progeny individuals with contrasted PgGs content (Figure 4), could be responsible for variation in anthocyanin biosynthesis.

The major colour-associated QTLs (PgGS, Anth, and COLOUR; ∼20% of explained variance) are located on the male M3A LG. The most likely underlying candidate gene is ANR. In the pseudo F₁ population studied, each individual is heterozygous at a given locus. Therefore, only one ANR homoeo-allele localised on the F. vesca–derived A subgenome underlies the male (“CF1116”) M3A colour-associated QTLs. Mapping the WGS of the “Capitola” and “CF1116” parents to the FvH4 F. vesca genome (Edger et al., 2018), which is syntenic to the subgenome Fvb3-4 (3A) (Hardigan et al., 2021a), allowed us to identify a large deletion (18 bp) in the 5′UTR of an ANR homoeo-allele (homoeo-allele g) carried by the male parent (“CF1116”). Combining genotyping by CE with visual scoring of fruit colour of the segregating population unequivocally associated intense red fruit colour to the presence of the ANR homoeo-allele g (Figure 6C).

These genetic indications, which support the contribution of ANR to the control of anthocyanin accumulation in strawberry, are in line with our previous molecular findings in F. vesca (Fischer et al., 2014). In the wild diploid F. vesca, ANR is encoded by a single gene localised on Fvb3. ANR converts the anthocyanidins pelargonidin and cyanidin to flavan-3-ols (epiafzelechin or epicatechin); concurrently, anthocyanidins can also be converted to anthocyanins (Figure 1A). In strawberry, flavan-3-ols accumulate in green fruit during the early stages of fruit development while the red-coloured anthocyanins accumulate after the onset of ripening (Fischer et al., 2014). Moreover, ANR is highly expressed in fruit cortex during early fruit development, up to the white fruit stage, and moderately so thereafter (Hawkins et al., 2017), suggesting a negative correlation between ANR transcript abundance and anthocyanin accumulation in the fruit. ANR may thus play a prominent role in the determination of fruit colour by controlling the trade-off between flavan-3-ols biosynthesis and anthocyanin biosynthesis. Indeed, we already demonstrated that silencing ANR in F. vesca fruit redirects phenypropanoid flux from flavan-3-ols to anthocyanins, resulting in early anthocyanin accumulation in the fruit (Fischer et al., 2014). Conversely, overexpression of tea (Camellia sinensis) ANR genes in tobacco results in a significant loss of flower red-pigmentation due to reduced ANTHc (Zhao et al., 2017). Altogether, these published results are consistent with our observation of a negative correlation between ANR transcript abundance and PgGs content in progeny individuals (Figure 4). They further support the hypothesis that the 18 bp ANR 5′UTR deletion found in the male homoeo-allele g may result in down-regulation of ANR, promotion of the flux from anthocyanidins to anthocyanins, accumulation of PgGS in the fruit and consequently more intense red colour in individuals of the progeny where homoeo-allele g is present.

The mechanism by which deletion in 5′UTR of ANR would achieve this effect remains to be elucidated. The 5′UTR plays a pivotal role in controlling gene expression through the regulation of transcript abundance or the alteration of mRNA translation efficiency or stability. Moreover, 5′UTR mutations can control transcript abundance by altering functional elements (Lim et al., 2021). The 5′UTR mutation in ANR homoeo-allele g deletes a cis-regulatory element which is a putative MADS box binding motif. By mining strawberry genomic databases, we found four MADS box genes that may possibly interact with the deleted motif and are expressed in fruit cortex along development (Hawkins et al., 2017). As the involvement of MADS-box in anthocyanin accumulation has been reported in different species including bilberry (Jaakola et al., 2010) and pear (Wang et al., 2017), which are also Rosaceae species, future work would focus on the possible involvement of the identified MADS box genes in the regulation of ANR and thus in the control of fruit colour in strawberry.