Deciphering the genetic control of fruit texture in apple by multiple family-based analysis and genome-wide association

Highlight A distinct set of QTLs related to mechanical and acoustic fruit texture features were identified in apple. Through a GWAS approach, the specific genetic control of these subtraits was elucidated.


Introduction
Fruit ripening is an orchestra of physiological changes occurring to render fruits more attractive and palatable (Giovannoni, 2001). This important quality feature depends on the dismantling of the primary cell wall polysaccharide complex by a series of cell wall-modifying proteins (Brummell and Harpster, 2001;Brummell, 2006). This synergic enzymatic action leads to different types of fruit texture in apple, from soft and mealy to firm and crisp, suggesting that rather than a single trait, fruit texture can therefore be considered as a multiple feature, with distinct specific mechanical and acoustic properties (Szczesniak, 2002;Zdunek et al., 2010;Costa et al., 2011). In the last decades, the control of fruit texture has represented a major goal towards the improvement of shelf-life performance (Matas et al., 2009). This aspect is of crucial importance, especially in the case of year-round fruit marketability and shipping overseas. The use of transgenic lines (Kramer and Redenbaugh, 1994) and whole-transcriptome platforms has in fact identified several genes involved in cell wall metabolism, such as those encoding expansin, pectin acetylesterase, xyloglucan endotransglycosylase, pectin methylesterase, pectate lyase, and polygalacturonase (Rose et al., 2002(Rose et al., , 2004Marín-Rodríguez et al., 2003;Eriksson, 2004;Brummell, 2006;Vicente et al., 2007;Janssen et al., 2008;Moore et al., 2008;Soglio et al., 2009;Costa et al., 2010a;Jimènez-Bermudéz et al., 2002;Segonne et al., 2014). In support of these studies, genomewide mapping of quantitative trait loci (QTLs) controlling fruit firmness is also compelling, with the final purpose of identifying the most valuable molecular markers suitable for marker-assisted breeding programs. An exhaustive knowledge of the fruit texture genetic make-up is in fact essential to guide the selection of the most valuable ideotypes in breeding by design approaches (Peleman and Van Der Voort, 2003). In this regard, it is worth emphasizing that the majority of QTLs associated with fruit texture characteristics have been focused mainly on one measurement, fruit firmness, and usually restricted to one or a few bi-parental mapping populations (Harada et al., 2000;King et al., 2000;Liebhard et al., 2003;Oraguzie et al., 2004;Kenis et al., 2008;Costa et al., 2010b;Chagné et al., 2014;Ben Sadok et al., 2015). To overcome this constraint, an important effort was represented by the simultaneous analysis of multiple populations connected in a common pedigree scheme, namely pedigreebased analysis (PBA) (Bink et al., 2008). This method has already been successfully employed to target QTLs in apple (Bink et al., 2014;Allard et al., 2016) as well as in cherry (Rosyara et al., 2013) and peach (Fresnedo-Ramírez et al., 2015. In addition to this, association mapping has also been widely employed as a complementary strategy to classical QTL mapping (Rafalski, 2010). Although this approach was initially employed in annual crops (Thornsberry et al., 2001;Weber et al., 2008;Stracke et al., 2009) and forest trees (Neale and Savolainen, 2004;González-Martínez et al., 2007Eckert et al., 2009;Neale and Kremer, 2011), it has also recently been exploited in fruit tree crops, such as grapevine (Cardoso et al., 2012), peach , and apple (Kumar et al., 2013(Kumar et al., , 2015. Especially in the latter species, a major QTL for fruit firmness was observed on chromosome 10, which coincided with MdPG1, a gene known to encode polygalacturonase playing a pivotal role in the depolymerization of pectins (Sitrit and Bennett, 1998;Brummell and Harpster, 2001). These investigations, however, are characterized by a low phenotyping resolution, to date recognized as the major operational bottleneck limiting the power of genetic analysis (Cobb et al., 2013).
To this end, the dissection of the fruit texture complexity was carried out with a texture analyzer, an instrument that has already demonstrated its reliability in dissecting the apple fruit texture into mechanical and acoustic signatures . To make advances in the deciphering of the genetic control of fruit texture in apple, a double approach was employed. Initially six full-sib pedigreed families were investigated through a PBA approach to detect QTLs associated with mechanical and acoustic fruit texture features. These regions were further complemented and validated by a genome-wide association study (GWAS) performed on a large apple germplasm collection to exploit a much larger range of both genetic and phenotypic variations.

Plant materials
In this study, two groups of plant materials were employed. The first was represented by a pedigree composed of 13 parental lines and six full-sib populations (for a total of 416 individuals). The scheme generated with PediMap (Voorrips et al., 2012;Supplementary Fig. S1; Supplementary Table S1 at JXB online) shows the maternal and paternal descendants from founders to progeny. The second group was represented by a collection of 387 apple accessions (Malus×domestica species; Supplementary Table S2), retrieved from the general apple variety repository available at the Fondazione Edmund Mach (FEM). The germplasm collection and two full-sib progeny (i.e. 'FjDe', 'Fuji'×'Delearly'; and 'FjPL', 'Fuji'×'Cripps Pink') were planted at the experimental orchard of FEM, while the other four full-sib families ('GDFj', 'Golden Delicious'×'Fuji'; 'GaPL', 'Gala'×'Cripps Pink'; 'GaPi', 'Gala'×'Pinova'; and 'FjPi', 'Fuji'×'Pinova') were chosen from the ongoing breeding program at the Laimburg Research Centre for Agriculture and Forestry. Within each group of plant materials (pedigreed full-sib families and germplasm), trees were at a full fruit-bearing stage at the time of the analysis. Moreover, while individuals from full-sib families were characterized by a single and original tree, the apple accessions included in the germplasm collection were represented by triplicates.
Fruit harvesting and high-resolution texture phenotyping Fruit, from both the pedigreed full-sib families and cultivar collection, were harvested at a commercial maturity stage according to typical fruit characteristics, such as starch degradation (selected at 7 on a scale of 1-10) and skin color. After harvest, fruit were stored for 2 months at 2 °C and 95% relative humidity. Prior to phenotyping, apples were maintained at 20 °C overnight.
The high-resolution phenotyping of fruit texture was carried out with a computer-controlled texture analyzer TAXTplus (StableMicroSystem, Godalming, UK), according to the protocol described by Costa et al. (2011Costa et al. ( , 2012. Since the texture analyzer was equipped with an AED (Acoustic Envelope Device), for each sample (flesh disc) a simultaneous profiling of the mechanical displacement and acoustic response was acquired. The combined profile was further processed with an ad hoc macro for the digital definition of 12 parameters, related to both mechanical (eight) and acoustic (four) texture properties (specified in Fig. 1). The phenotype was assessed for two consecutive experimental years and represented by BLUP (Best Linear Unbiased Prediction, 'Ime4' R package), which adjusted the mean by reducing the error variance.

SNP genotyping
Genomic DNA was isolated from young leaves with the Qiagen DNeasy Plant Kit (Qiagen). DNA quantity and quality were measured with a Nanodrop ND-8000 (ThermoScientific, USA). The genotyping was carried out with the 20K Infinium array (Illumina; Bianco et al., 2014;www.fruitbreedomics.com). The SNP data, initially processed with the GenomeStudio Data Analysis software, were finally elaborated with ASSIsT (Di Guardo et al., 2015), a dedicated stand-alone pipeline to filter and re-edit SNP calls with a distorted segregation pattern due to the presence of null alleles (Pikunova et al., 2013).

Linkage mapping and SNP marker consensus order
The array of SNPs segregating within the six bi-parental pedigreed populations was used to create a bi-parental integrated map for each full-sib family, using the software JoinMap 4.1 (Van Ooijen, 2006). Markers were initially clustered in linkage groups with a minimum LOD value of 3.0 and further ordered with a recombination frequency of 0.45 and Haldane mapping function. The six genetic maps were finally merged into a consensus map using BioMarcator software v4.2 (Sosnowski et al., 2012), through the implementation of the ConsMap module. The new reference and harmonized marker order (Supplementary Table  S3) was exploited as a map input file for both PBA and GWAS analysis.

QTL mapping by pedigree-based analysis
The identification and mapping of QTLs at the genome-wide scale was carried out with FlexQTL™ (Bink and van Eeuwijk, 2009;Bink et al., 2014;www.wur.nl/en/show/flexqtl.htm). The Bayesian approach on which the software is based compares different models, considered as a random variable. The linear model is expressed as follows: where y is the observed phenotype, μ is the phenotypic mean, W is a matrix of a vector of regression on the QTL covariates (a), and e is the residual error of the model. The model operated by FlexQTL™ is based on independent assignment of QTL alleles to founders; therefore, the genotype of several individuals is a priori unknown. Since the joint posterior distribution for the number of QTLs cannot be analytically computed due to the high number of genotype combinations, a Markov chain Monte Carlo (MCMC) simulation is used. The effective chain size (ECS; Sorensen and Gianola, 2002), used to assess the sensitivity of the posterior inference, was considered statistically significant with values >100 for phenotypic mean (μ), QTL explained variance (σ 2 a ), QTL residual variance (σ 2 e ), and number of expected QTLs on the a priori distribution (N QTL ), respectively. In the analysis carried out here, 500 000 iterations were performed and a thinning of 500 was applied to reduce computation storage. For each run, the convergence quality was represented by a trace plot for visual inspection ( Supplementary Fig. S2). For each chromosome, the number of QTL(s) was inferred comparing models In the loading panels, red-colored arrows indicate the mechanical parameters, while acoustic parameters are shown with blue arrows. Each parameter is coded with a number as follows: 1, initial force; 2, maximum force; 3, final force; 4, mean force; 5, area; 6, force linear distance; 7, Young's modulus; 8, number of force peaks; 9, maximum acoustic pressure; 10, mean acoustic pressure; 11, acoustic linear distance; 12, number of acoustic peaks.
with an increased number of QTLs. In FlexQTL™, the most likely number of QTLs was inferred based on Bayes factors (BFs), which represent the ratio of the marginal likelihood under one model compared with the marginal likelihood under a second model. A 2ln(BF) ≥2, 5, or 10 indicates positive, strong, or decisive evidence of a model, respectively. Moreover, each QTL is considered as biallelic, with three possible genotypic conformations: 'QQ', 'Qq', and 'qq'. The analysis was carried out with an additive genetic model, assigning to 'QQ' and 'qq' a value of 1 and -1, respectively, while 'Qq' was equal to 0. The QTL genotype of each individual included in the pedigree is a priori unknown and the alleles are assigned to founders tracing their transmission to offspring. To reduce computation time and increase marker informativeness, the initial set of 10 695 markers was finally converted into 1045 haploblocks with PediHaplotyper (Voorrips et al., 2016) Linkage disequilibrium analysis The SNP markers exploited on the germplasm collection were initially employed to estimate the linkage disequilibrium (LD) decay. From the 20K SNPs, 7378 markers were actually used for LD analysis, excluding rare alleles with a minor allele frequency (MAF) <0.05 and those showing an incongruent physical/genetic position. The pair-wise r 2 between SNP markers was calculated with Plink (Purcell et al., 2007). For both a chromosome-wise and genome-wide scale, the LD decay was depicted by plotting the pair-wise r 2 value against the corresponding physical distance on the genome (bp). The estimation of the LD decay distance was defined by crossing the r 2 baseline (based on the 95th percentile of the marker distribution, according to Breseghello and Sorells, 2006) and the locally weighted polynomial regression-based fitting curve (LOESS) fitted to the plot ('stats' R package). For each chromosome, the LD level was also depicted by partitioning the chromosomal regions into segments of strong LD with Haploview (Barrett et al., 2005).

Population structure and genome-wide association mapping
The level of genetic stratification was assessed with STRUCTURE v2.3.1 (Pritchard et al., 2000). To this end, 17 simple sequence repeats (SSRs; Supplementary Table S4) were amplified according to the protocol reported in Di Guardo et al. (2013). The SSR genetic data were further used to compute the posterior probability [Pr(X|K)], given a specific number of group K (ranging from K=2 to K=8). The computation was carried out performing five independent runs of 1 000 000 burn-in generations and considering the admixture model. The most probable number of populations was identified with STRUCTURE HARVESTER (Earl and vonHoldt, 2012), and the final population structure matrix (Q) was further implemented as a covariate in GWAS analysis.
The marker-trait association analysis was carried out with TASSEL v3 software (Bradbury et al., 2007; http://www.maizegenetics.net). The significance of the association was tested implementing both the general linear model (GLM) and the mixed linear model (MLM). The GLM (Pritchard et al., 2000) was computed correcting for population structure. The MLM (Zhang et al. 2010), instead, included both fixed and random effects, allowing the incorporation of genetic relationship as follows: where y represents the phenotype (vector of observation), β is an unknown vector containing fixed effects (marker and Q), u is an unknown vector of random additive genetic effects, X and Z are the known designed matrices, and e is the vector of random residuals. The MLM considered both Q (population structure) and K (kinship) matrixes as covariates for population and parental relationship correction (false positive). Significant associations were selected according to a P-value ≤0.05, after false discovery rate (FDR) correction for multiple comparison according to the procedure of Benjamini and Hochberg (1995) using the 'stats' R package. For each trait considered in the association, the choice of the model was suggested by the visual inspection of the Q-Q plot, obtained with the 'qqman' R package.

Gene expression analysis by RT-qPCR
To assess whether a change in gene expression corresponds to a different QTL estimated genotype, the transcription profiling of MdACO1 and MdPG1 was assessed at harvest and after 2 months of cold storage. To achieve this goal, the three parental lines ('Delearly', 'Fuji', and 'Cripps Pink') together with four seedlings for 'FjDe' (45, 125, 10, and 14) and 'FjPL' (23,25,35,and 68) were selected. Fruit mesocarp was cut, frozen in liquid nitrogen, ground into a fine powder, and stored at -80°C until processing. RNA extraction, quantification, and reverse transcription-quantitative PCR (RT-qPCR) were carried out according to the methods described in Busatto et al. (2015). The final Ct is represented by the average of three independent normalized expression values for each sample, and an actin gene (MdACT) was employed as housekeeping gene (Di Guardo et al., 2013). For each gene, a pair of discriminant and specific primers was designed (Supplementary Table S5), using Primer3 (http://primer3.ut.ee) and Primique (http://cgi-www.daimi.au.dk/cgi-chili/primique/front.py).

High-resolution phenotyping of fruit texture behavior in apple
The apple fruit texture was assessed with a novel and sophisticated texture analyzer (Costa et al., 2011. The overall phenotypic fruit texture variability was initially represented by a principal component analysis (PCA) plot ( Fig. 1). For both groups of plant materials, the fruit texture parameters were similarly oriented, with a consistent incidence of the two principal components (PCs) chosen to define the PC hyperspace. Comparing the two plots, PC1 explained 71.6% of the total phenotypic variance in the pedigreed full-sib families (Fig. 1A, B) and 79.6% in the germplasm collection (Fig. 1C, D), while PC2 accounted for 19.9% and 12.7%, respectively. Individuals were distributed following the loadings' projection (Fig. 1B, D) represented by the 12 texture parameters, which clearly discriminated the two signatures. In both scenarios, the variables were in fact distinctly oriented towards two PCA quadrants. The mechanical parameters (highlighted by the numerical code from 1 to 8) were mostly plotted in the PC1 positive/PC2 negative quadrant, while the acoustic parameters (9-12) were projected on the quadrant defined by positive values for both PCs, besides the number of force peaks (8). Although this index is considered a mechanical parameter, it is correlated more with the group acoustic indices, justified by the mechanism behind the generation of the acoustic response and pressure progression (Vincent, 1998).

Fruit texture QTL discovery through PBA
Each parameter obtained from the phenotypic dissection of fruit texture was finally exploited in marker-trait association studies. In the attempt to map the QTLs related to this feature, the Bayesian approach was initially employed. QTLs were identified and mapped on 13 chromosomes, on which the posterior QTL intensity exceeded the posterior probability threshold [2ln(BF) >2; Table 1]. The overall genome-wide Table 1

. For each trait assessed, the variance (Var), the mean (Mean), and the linkage group (LG) on which the QTL is identified are reported
For each QTL, the interval, the length, the mode (in cM), the 2ln Bayes factor (BF) for the presence of one or two QTLs (1/0 and 2/1, respectively), the probability (Prob), the additive effect size (AEt1), the additive variance (AVt1), and the weighted additive variance (wAVt1) are also shown. Each parameter refers to a haploblock, specified in the last column  Fig. S3) distinguished specific probability profiles for the two groups of texture-related parameters, acoustic and mechanical. For simplicity, the QTL differences are highlighted in Fig. 2, comparing the profiles of the maximum force (mechanical) with the number of acoustic peaks (acoustic). Chromosomes 3 and 10 showed QTLs commonly shared by both features (Fig. 2; Supplementary  Fig. S3). For the maximum force, in particular, the QTL mapped on chromosome 3 ( Fig. 2A) is located on a single genomic region [2ln(BF) 1/0 =13.6], with a mode at 55 cM and an allelic effect (AEt1) of 1.54 (Table 1). This position was also similar to the rest of the mechanical parameters, spanning from 55 cM to 57 cM (Supplementary Fig. S3; Table  1). The only difference was observed for the number of force peaks, which showed the QTL at 10 cM. This observation, however, additionally confirms the association of this parameter with the group of acoustic parameters. In the case of the acoustic peaks (Fig. 2B), two QTLs were instead suggested [2ln(BF) 2/1 =2.7], with a mode at 3 cM and 28 cM and an AEt1 of 10.55 and 9.39, respectively. Among the acoustic subtraits, the QTL on chromosome 3 was also identified for the acoustic linear distance, but at 28 cM. On chromosome 10, a single QTL associated with the number of acoustic peaks (Fig. 2B) was shown [2ln(BF) 1/0 =10.6] and located at 40 cM with an AEt1 of 18.25. This position was also similar to the rest of the acoustic parameters (at 40 cM and 42 cM), beside the mean acoustic pressure that showed the QTL peak at 35 cM, as was also observed for the number of force peaks (Table 1). For the maximum force, two QTLs were instead observed on this chromosome [2ln(BF) 2/1 =4.2]. The first was located at 20 cM, with an AEt1 of 0.91, while the second was mapped at 45 cM with an AEt1 of 1.92. These two regions were also consistent across the mechanical parameters (spanning between 19 cM and 20 cM for the first QTL and 44-46 cM for the second), with two exceptions. The force linear distance in fact showed only one QTL (with a low probability and effect) at 49 cM, while for the Young's modulus (or elastic modulus) no QTL was observed (Table 1; Supplementary Fig. S3). Beside these, other QTLs were identified with a more specific pattern. The two major genomic intervals showing QTLs associated with mechanical parameters were mapped on chromosomes 11 and 16 ( Fig. 2; Supplementary Fig. S3). The QTLs positioned on chromosome 11, across the several mechanical parameters, were located between 41 cM and 49 cM (Table 1). On this chromosome it is also interesting to note the QTL positioned at 14 cM [2ln(BF)=7.9 and AEt1 of 0.11] and related to the Young's modulus. The different and original positioning of this QTL within the class of mechanical parameters can be due to the fact that the Young's modulus depends on the elasticity of the sample (the ratio between stress and strain) rather than fruit firmness. The particular behavior observed for the Young's modulus, besides its projection over the PCA plot (Fig. 1B, D), is moreover validated by the QTL profile detected on chromosome 16. As reported for chromosome 11, this QTL is also specifically associated with the mechanical parameters ( Supplementary Fig. S3), with a mode located at ~32 cM, with the exception of the Young's modulus, where this QTL was not detected (similarly to chromosome 10). In the second subset of QTLs (related to the acoustic parameters), several genomic regions located on chromosome 1 and associated with mean and maximum acoustic pressure, number of force peaks, as well as the Young's modulus were identified. Regarding force peaks, other QTLs were moreover mapped on chromosomes 5, 8, and 9. This latter chromosome, together with 13 and 14, was also associated with the number of acoustic peaks ( Fig. 2B; Table 1; Supplementary Fig. S3). The simultaneous presence of QTLs detected on chromosome 9 for both number of force and acoustic peaks strengthens the relationship between these two parameters. Especially for the number of acoustic peaks, chromosome 9 and 14 were found to be the most important, showing an allelic effect of 10.27 and 16.48, respectively (Table 1). In particular, the QTL on chromosome 14 is characterized by an estimated genotype (Fig. 3A) consistent with the acoustic performance (assessed as the number of acoustic peaks) of the six parental cultivars (Fig. 3B, C). Among the group, 'Pinova' and 'Fuji' were distinguished by the highest acoustic response, as depicted in the 2D-PCA plot (Fig. 3B) and loading projection (Fig. 3C). The superior crispness performance of these two apple cultivars is also confirmed by the homozygous state of the positive estimated QTL allele ('QQ'). In contrast, cultivars with a mealy texture, such as 'Delearly' and 'Royal Gala', are plotted on the other extreme of the 2D-PCA plot, showing a 'qq' genotype for this QTL. The effect of the estimated allele for the QTL on chromosome 14 was further investigated on the six progeny (Fig. 4). FlexQTL™ estimated a 'QQ' genotype for 'Fuji' and 'Pinova', a heterozygous 'Qq' genotype only for 'Golden Delicious', and a 'qq' genotype for the other three varieties ('Delearly', 'Royal Gala', and 'Cripps Pink'). Taking into account that the 'Q' allele is known to increase the phenotypic performance, it is worth noting that only the seedlings of 'FjPi' and half of the progeny of  'GDFj' (with a 'QQ' genotype) were distinguished by the highest acoustic response (Fig. 4), underlying the role of this QTL in the control of the acoustic properties in apple.

Analysis of the linkage disequilibrium in domesticated apple
The level of LD was determined to verify the genetic associations between loci and to scan the LD decay over each chromosome. To this end, from the total set of 10 695 selected SNP markers, 7378 were effectively employed in the computation, after subsequent filtering steps. From the initial data set, besides SNPs with a MAF <0.05, markers showing an inconsistent position between the physical location on the genome and the consensus genetic map were also excluded ( Supplementary Fig. S4). In the germplasm collection investigated here (and represented by 387 apple accessions), the LD decay was estimated to extend for an average up to r 2 =0.19 at the genome-wide level (Fig. 5A) corresponding to ~400 kb, and spanning from a maximum of r 2 =0.28 for chromosome 16 (Fig. 5B) to a minimum of r 2 =0.13 for chromosome 17 (Fig.  5C). In parallel, the presence of distinct LD blocks over the genome was illustrated with an LD heatmap ( Supplementary  Fig. S5), highlighting specific genetic fixation for each apple chromosome. Of all chromosomes, chromosome 16 is characterized by the highest LD value, showing an LD block of 2675 kb (Fig. 5B). These results indicate that the extent of LD in apple is shorter than in peach  but larger than that in other species (e.g. grapevine; Myles et al., 2009).

Fruit texture genetic dissection by GWAS
The QTLs identified with the PBA approach were further complemented by GWAS. From the entire germplasm collection assessed to estimate the LD decay, 233 accessions were used for both population structure and marker-trait association. Individuals were assigned to three subpopulations (K=3; Supplementary Fig. S6) following the plateau criterion (Falush et al., 2007), the non-parametric Wilcoxon test (Rosenberg et al., 2002), and the ΔK method proposed by Evanno et al. (2005). Beyond this point, the mean log-likelihood values tend to a plateau together with an increased SD, which became clearly evident from K=5.
From the 20K SNPs present in the array, 10 558 were finally exploited in the GWAS computation, performed with the MLM implemented in TASSEL v3.0. As a first attempt to dissect the genetic control of apple fruit texture, the two principal components (PC1 and PC2) were initially implemented as phenotypic traits. In this case, PC1 was considered to capture the overall texture variability, the entire group of parameters being oriented towards its projection (Fig. 1C, D). PC2, instead, was employed to discriminate the mechanical from the acoustic subset of variables. The MLM module identified for PC1 a major QTL on chromosome 10 (Supplementary   Supplementary Fig. S7A) For each SNP, the chromosome on which the marker is mapped (Chr), the genetic position (cM), the P-value, the name (SNP), and the percentage of phenotypic variability explained (r 2 ) are indicated  Fig. S7A) with 10 markers exceeding the FDR-corrected P-value threshold. The phenotypic variability explained by the markers spanned from 9% to 13% and were found positioned on the consensus map between 42 cM and 47 cM (Table 2). Within this marker set, it is worth highlighting five SNPs, coded as FEM_cg_9,11,17,18,and 19, which are custom SNPs specifically designed on polymorphisms discovered on re-sequencing the full length of MdPG1, a gene playing a key role in the fruit softening process in apple (Wakasa et al., 2006;Costa et al., 2010b;Longhi et al., 2012). In particular, FEM_cg_19, also named Md-PG1SNP (Baumgartner et al., 2016), is an SNP highly correlated with the microsatellite marker Md-PG1 SSR 10kd, previously associated with the fruit texture behavior in apple (Longhi et al., 2013). In contrast, when PC2 was used in the association analysis not a single SNP was identified as statistically significantly associated ( Supplementary Fig. S7B). It is moreover worth noting that when single texture subtraits were used as phenotype, the association result was consistent with the profile obtained for PC1, as shown in Supplementary Fig. S7C and D for the maximum force and number of acoustic peaks, respectively. Although the two groups of variables are oriented towards two different PCA quadrants, they are, however, commonly projected along the PC1 orientation (Fig. 1D). Within the panel of apple accessions employed here, PC1 accounts for 79.6% of the total phenotypic variance, thus influencing the genetic association of each single texture parameter. PC2, instead, is orthogonally oriented with respect to PC1 and more related to the difference between the two groups of variables, and therefore more effective in the dissection of the genetic control of the two texture signatures.
To decipher more specifically the genetic regulation of fruit texture properties, a second round of association was performed. Since it is already well known that crisp apples are more appreciated by consumers, from the initial set of accessions used in GWAS, genotypes distinguished by the unfavorable homozygous allelic configuration for MdPG1 and MdACO1 were removed. The effect of these two genes on the fruit texture in apple largely depends on the interaction of the physiological processes they control. MdPG1 is involved in the dismantling of the cell wall/middle lamella structure (Brummell and Harpster, 2001;Brummell, 2006) and its effect in apple seems to be more relevant than in other climacteric species. In tomato, in fact, the role of this gene alone does not impact the fruit texture physiology significantly (Sheehy et al., 1988;Smith et al., 1988;Giovannoni et al., 1989). The other gene, MdACO1, regulates the last step of ethylene biosynthesis (Bleecker and Kende, 2000). Although in climacteric fruit the amount of this hormone is known to control several processes (Rose et al., 1998;Costa et al., 2005;Wakasa et al., 2006), the co-existence of ethylene-dependent and -independent regulation of fruit texture also has been proposed, as demonstrated in melon (Nishiyama et al., 2007) as well as in apple (Tadiello et al., 2016). This dual mechanisms can explain why QTLs for fruit firmness have been collocated with MdPG1 but not MdACO1. To investigate the consistency between the QTL genotypes estimated by FlexQTL™ and the expression of a gene included in the corresponding genomic interval, the transcript accumulation of MdPG1, together with MdACO1, was assessed. The transcript profiling was carried out in two groups of seedlings chosen between 'FjDe' and 'FjPL' populations. Among the four individuals selected in the first population, two (FjDe_10 and FjDe_14) were characterized by a 'QQ' genotype estimated for the QTL on chromosome 10 (and coincident with the genetic position of MdPG1), while the other two (FjDe_45 and FjDe_125) were distinguished by a 'Qq' genotype. These two QTL genotypes are, moreover, consistent with the allelotype configuration of Md-PG1SNP. The 'Q' and 'q' alleles are in fact linked to the 'C' and 'T' allele of this marker, in agreement with the genotype of the parental lines. 'Fuji' ('QQ' estimated genotype) is indeed distinguished by a 'CC' allelic state for Md-PG1SNP, while 'Delearly' has an 'TC' allelotype. The 'T' allele (related to the 'q' QTL allele) therefore segregates within the 'FjDe' population, contributing to a loss of fruit firmness. This association is further validated by the Pearson correlation value (R 2 -0.8) between MdPG1 expression and the fruit firmness assessment depicted in Fig. 6. In both groups of parental lines (group 1 in Fig. 6) and 'FjDe' individuals (group 2), it is clear that high loss of firmness corresponds to high MdPG1 expression. The genotypes distinguished by a 'q' allele (FjDe_45, FjDe_125, and the parental variety 'Delearly') are in fact characterized by a considerable MdPG1 expression already at harvest. This observation was additionally confirmed by the analysis carried out on the second population, 'FjPL'. Since the two parents ('Fuji' and 'Cripps Pink') do not segregate for this QTL, all four seedlings (group 3) are characterized by a 'QQ' estimated genotype (MdPG1SNP_CC). Due to the role of these two genes, the breeding activities oriented towards the selection of firm and crisp apples no longer consider cultivars with these unfavorable genotypes as parental lines. Although this second panel of accessions is composed of only 64 individuals, it captures the real phenotypic variability used nowadays by breeders. This second GWAS was carried out with a GLM, selected on the basis of the Q-Q plot inspection. The result of this re-shaped phenotypic variance, obtained by fixing the effect of the two loci, was evident in the association analysis depicted in Fig. 7 and computed for the maximum force (Fig. 7A) as well as the number of acoustic peaks (Fig. 7B). In both associations, no major QTL on chromosome 10 was detected, especially in the case of the mechanical parameter. In contrast, when the number of acoustic peaks was considered, other regions were identified and located on chromosomes 2, 14, and 15 ( Fig. 7B; Supplementary Table S6). These results provide evidence about the distinct genetic control for the two texture properties in apple, and suggest the role of chromosome 14 in the determination of acoustic properties, as underlined by the PBA results.

QTL anchoring and in silico gene annotation
To investigate further the role of the QTL mapped on chromosome 14 in the regulation of the acoustic component of fruit texture, an in silico search and annotation of the candidate genes included in the genomic interval, together with their transcription profiling, was carried out. Genes were searched and annotated within an interval of 400 kb (LD block) up and downstream from the QTL peak determined by both the PBA and GWAS approach (Table 3). Among them, it is important to highlight categories of importance. The first is represented by proline-rich proteins (PRPs). This class of cell wall-modifying proteins (CWMPs) incluses proline and hydroxyproline peptides and seems to be involved in the cell wall metabolism of several species, such as cotton, carrot, and Arabidopsis (John and Keller, 1995;Fowler et al., 1999;Holk et al., 2002). It is also interesting to note that PRP-related genes are expressed in immature watermelon fruit (Guo et al., 2011), and therefore not related to the late ripening dismantling process of the cell wall. Another important CWMP class is represented by expansin, a type of protein involved in the architectural re-modeling of the cell wall causing a disruption of the non-covalent bonds between the hemicellulose matrix and cellulose microfibrils (Cosgrove, 2000), exposing the cell wall structural polymers to the action of other CWMPs. Within the QTL computed with GWAS, a xylanase (xylo-glucangalactosyltransferase) was also found. Heteroxylans are a divergent group of polymers contributing to the cell wall structure, although in dicots they are less abundant than xyloglucan (Johnston et al. 2013).

Table 3. Gene annotation within the QTL interval mapped on chromosome 14 and identified through both GWAS (upper part) and PBA (lower part) approaches
For each gene, the NCBI and GDR gene ID, gene function, and conting co-ordinates are reported. In addition, the closest SNP with its relative genomic information is also provided Gene ID Xylanases are thus involved in the cellulose-non-cellulosic framework, and their role has already been reported in fruit. This gene was in fact expressed during the fruit softening process in papaya, according to an ethylenedependent pattern (Manenoi and Paull, 2007). Finally, it is also worth mentioning the identification of a fucosyltransferase, an enzyme involved in xyloglucan metabolism, as documented in Arabidopsis (Rocha et al., 2016).
To assess further the mode of action of these genes in the regulation of fruit texture, the whole transcriptomic survey presented by Tadiello et al. (2016) was re-examined. In this particular case, the expression profile of two expansins (MDP0000423907 and MDP0000193025) and a fucosyltransferase (MDP0000230681) was selected from the analysis carried out with the whole-genome apple array (WGAA). The expression analysis of these elements was retrieved from three specific samples of the reference apple cultivar 'Golden Delicious', at harvest and after 1 week of post-harvest shelflife ripening under both normal and 1-methylcyclopropene (1-MCP)-treated conditions ( Supplementary Fig. S8). As described in Tadiello et al. (2016), the control post-harvest sample is characterized by an important loss of acoustic performance with regards to harvest, while the application of the ethylene competitor (1-MCP) effectively limited the cell wall degradation. Interestingly, the genes identified here and associated with the regulation of the acoustic component of texture are not stimulated by ethylene and do not participate in the major cell wall-degrading events, since their expression does not change from harvest to post-harvest. One expansin, in particular (MDP0000423907), is induced when the acoustic performance is promoted (1-MCP treatment), meaning that its role, rather than in the dismantling process leading to softening, is an involvement in the maintenance of the cell wall architectural structure.

Conclusion
Fruit texture in apple is made up of multiple subtraits, most of which (especially the acoustic ones) are poorly investigated, in particular for genetic purposes. So far, this limitation led to the identification of markers suitable to assist in the selection of fruit firmness only, although, for apple, the feature most preferred by consumers is crispness. The coupling of PBA with GWAS enabled the genetic deciphering of fruit texture control, identifying important QTLs associated with both texture features. The comparison of the results obtained by the two genetic approaches highlighted an inventory of genomic intervals specifically associated with mechanical and acoustic parameters, respectively, hypothesizing that these subtraits are effectively controlled by different genetic mechanisms. In the near future, in the new high quality breeding materials, the alleles of the markers currently in use will be quickly fixed as a result of recurrent rounds of marker-informed parental selection (MAPS; marker-assisted parent selection) and the subsequently assisted selection of seedlings (MASS; marker-assisted seedling selection). The information presented here can therefore be taken into consideration to design novel markers useful to identify novel apple accessions distinguished by superior fruit crispness.

Supplementary data
Supplementary data are available at JXB online. Fig. S1. Pedigree structure of the six full-sib families. Fig. S2. QTL trace plot position output of FlexQTL™ showing the convergence of each single run. Fig. S3. Overall QTL probability profile defined for each parameter over the entire genome. Fig. S4. For each chromosome, the correlation between the physical and genetic position of each SNP marker employed in the analysis is shown. Fig. S5. LD decay plot and Haploview LD pattern for each apple chromosome. Fig. S6. Inferred population structure. Fig. S7. Genome-wide association scan of the association between SNP marker loci and PC1 (A), PC2 (B), maximum force (C), and number of acoustic peaks (D) computed with the collection of 233 apple accessions. Fig. S8. Expression profile of three genes (MDP0000423907 in blue, MDP0000193025 in red, and MDP0000230681 in green) in three samples of 'Golden Delicious', according to the work of Tadiello et al. (2016). Table S1. Mating scheme of the six full-sib pedigreed families Table S2. List of the apple accessions employed in the LD analysis. Table S3. List of SNP markers employed in this study. Table S4. List of 17 SSR markers used to analyze the population structure. Table S5. List of primer for RT-qPCR analysis. Table S6. List of SNPs identified in the GWAS analysis carried out with the pre-selected group of accessions.