Association genetics of chilling injury susceptibility in peach (Prunus persica (L.) Batsch) across multiple years

Peach and nectarine (Prunus persica L.) are highly perishable; they ripen and deteriorate quickly at ambient temperature. Storage at low temperature (0–5 °C) is a common strategy used to slow the ripening processes and to extend shelf life. However, if susceptible varieties are held too long at a low temperature, they will not ripen properly and will develop chilling injury (CI) symptoms like mealiness (M), flesh browning (FB), and flesh bleeding (FBL). Understanding the genetic control of these traits to produce CI resistant cultivars will greatly benefit producers, shippers and consumers. In this study, we evaluated a population of 51 individuals from Pop-DG across 4 years with CI traits observed in one or two time points to detect molecular marker association with selected 960 single-nucleotide polymorphisms (SNPs) from 1,536 SNPs chip. Genotypic and phenotypic data were analyzed by general linear model and mixed linear model to see comparative results from both analyses. Among 960 SNPs used, 22 SNPs were found associated with CI susceptibility traits like M, FB, and FBL. Many SNP markers were located in or close to previously reported quantitative trait loci mapped by linkage analysis.


Introduction
Peach tree (Prunus persica) is a species of Prunus, a genus that also includes nectarine, plum, apricot, cherry, and almond belonging to the subfamily Prunoideae of the family Rosaceae. It is considered one of the genetically most well characterized species in the Rosaceae, and it has distinct advantages that make it suitable as a model genome species for Prunus as well as for other species in the Rosaceae (Abbott et al. 2002;Shulaev et al. 2008). Peach is a diploid with n = 8 and has a comparatively small genome currently estimated to be *220-230 Mbp based upon the peach v1.0 assembly. Peach has a relatively short juvenility period of 2-3 years compared with most other fruit tree species that require 6-10 years. Peach and nectarine (Prunus persica L.) are highly perishable; they ripen and deteriorate quickly at ambient temperature. Storage at low temperature (0-5°C) is a common strategy used to slow the ripening processes and to extend shelf life. Susceptibility of stone fruit to chilling injury (CI) is highly influenced by the genetic background of the cultivar (Peace et al. 2006). The physiological basis of CI symptoms has been studied in detail and well reviewed in peach (Lurie and Crisosto 2005). However, the exact mechanism by which CI affects a commodity is not fully understood. It has been shown to involve loss of membrane integrity and ion leakage from cells and changes in enzyme activity (Brummell et al. 2004a, b), but exactly why some crops are susceptible and some resistant still remains unclear.
Understanding the genetic control of these traits, to grow only cultivars free of CI susceptibility, promises to greatly benefit producers, shippers, and consumers in the peach industry. The major symptoms of CI are flesh mealiness (M), flesh browning (FB), and flesh bleeding (FBL) (Crisosto et al. 1999). M is a fruit flesh textural disorder, where affected ripe fruit have a dry, grainy feel when chewed. In simple terms, mealy fruit are dry and soft when ripe (Ju et al. 2000). FB is often seen in mealy fruit, although it can also occur in the absence of mealiness (Crisosto et al. 1999). It occurs when enzymes such as polyphenol oxidase act on phenolic substrates with which they come in contact. FBL is caused by the movement of water-soluble red pigments, probably anthocyanins, through the fruit flesh during cold storage or after subsequent ripening (Lurie and Crisosto 2005).
In the last decade, several linkage maps obtained using molecular markers have been constructed for peach (Aranzana et al. 2002;Dirlewanger et al. 2006;Etienne et al. 2002). A consensus map from an inter-specific almond 9 peach F2 population ('Texas' 9 'Earlygold', T 9 E) is considered the reference map of the Prunus genus (Howad et al. 2005). However, many important agronomic characters of Prunus species have not yet been mapped, and very few of those already mapped (such as major genes for disease and pest resistances, self-incompatibility, and several fruit quality traits such as flesh color, endocarp staining, flesh adherence to stone, non-acid fruit, skin pubescence, skin color, and fruit shape) are currently being used for marker-assisted selection (MAS) (Dirlewanger et al. 2004). The application of next generation sequencing technologies and bioinformatic scripts to generate high frequency single nucleotide polymorphisms (SNPs) distributed throughout the peach genome for use in genome mapping and phenotype selection and development of high density genetic linkage maps using SNP markers were considered for two breeding populations, Pop-DF ('Dr. Davis' 9 'F8, 1-42') with 117 progeny and Pop-DG ('Dr. Davis' 9 'Georgia Belle') with 64 progeny (Martínez-García et al. 2012). New, integrated, and saturated SNP linkage map and high density QTL discovery in our previous study and marker-trait association across 4 years in our present one will be a valuable resource for researchers working in peach and other related species. Transcriptomic analyses of two peach cultivars namely Oded and Hermoza, which differ in their resistance to CI, were examined after 2 weeks of cold storage at 5°C using ChillPeach cDNA microarray platform and identified 107 CGs proposed to be involved in CI (Dagar et al. 2012).
Genome-wide association analysis (GWAS) is a powerful approach to identify the causal genetic polymorphisms underlying complex traits (Riedelsheimer et al. 2012;Zhao et al. 2011). Development of genomics has provided alternative tools to improve breeding efficiency in plant breeding programs. Molecular markers linked to the causal genes and/or quantitative trait loci (QTLs) can be used for MAS (Xu and Crouch 2008). The key constraint for the successful use of association analysis in plants is the population structure and genetic relatedness, which can result in spurious marker-trait associations that may make it difficult to distinguish loci that truly affect the target traits (Chan et al. 2011;Ersoz et al. 2007;Gupta et al. 2005). Recent advances in genome sequencing and SNP genotyping have increased the applicability of association analysis for QTL mapping in crops (Morgante and Salamini 2003;Rafalski 2010).
SNPs have a low mutation rate and are evolutionarily stable from generation to generation across the genome (Batley and Edwards 2009). SNP markers have several advantages for genetic mapping over other molecular markers. SNPs have fewer detection/evaluation errors than simple sequence repeats (SSRs) (Hamblin et al. 2007) and map QTLs with greater precision than is possible with restriction fragment length polymorphisms (RFLPs) or SSRs (Yu et al. 2011). SNPs are often transferable across species within a genus (Grattapaglia et al. 2011). The development of a 6,654 peach SNP panel from which a 1,536 SNP set was selected for our previous study (Ahmad et al. 2011). Single nucleotide polymorphism (SNP) markers covering the entire genome are needed to enable molecular breeding efforts such as genome-wide association studies, fine mapping, genomic selection, and MAS in peach [Prunus persica (L.) Batsch] and related Prunus species. To date only a limited number of genetic markers, including SSRs have been available, this issue was additionally addressed by an international consortium (The International Peach SNP Consortium; IPSC) in a coordinated effort to perform genome-scale SNP discovery in peach using next generation sequencing platforms to develop and characterize a highthroughput Illumina Infinium Ò SNP genotyping array platform (Verde et al. 2012).
In previous studies, both major and minor QTLs controlling mealiness, browning, and bleeding were localized, using phenotypic data collected for three harvest seasons. The endopolygalacturonase gene, at the F-M locus, is responsible for a major QTL controlling both mealiness and bleeding, while one of the minor QTLs for bleeding was located close to the flesh color (Ogundiwin et al. 2007). QTLs for mealiness, graininess, leatheriness, and bleeding were localized on the LG4 of V 9 BT population and a major QTL for mealiness was validated in V 9 BT population and QTLs for browning were not found on LG4 in our previous work (Cantín et al. 2010). Our previous study resulted in first genetic linkage map of CI susceptibility in peach with SSR and SNP markers (Dhanapal et al. 2012). The aim of this study was to identify SNPs association with CI susceptibility in peach across multiple years. This study will represent the most complete and comprehensive association analysis to date of CI symptoms in Pop-DG progeny populations and contributes to our present understanding of the genetic control of CI symptoms in peach. Strong marker-trait association was detected for M, FB, and FBL making breeding for resistance or low susceptibility to CI an achievable goal. New, integrated, and saturated SNP linkage map and high density QTL discovery in our previous study and marker-trait association across 4 years in our present one will permit development of new markers to aid selection of new cultivars with low susceptibility or resistance to CI. To our knowledge this is the first GWAS in peach using SNP markers for CI across 4 years with one or two time points.

Materials and methods
Plant materials, DNA isolation, and quantification 'Pop-DG', a peach intraspecific cross between 'Dr. Davis' (female parent) and 'Georgia Belle' (pollen parent), was used in this study. 'Dr.Davis' is a modern canning peach cultivar and 'Georgia Belle' is a century-old fresh market peach cultivar which contrasts for many fruit quality and other CI related traits. 'Pop-DG' was created and managed at Kearney Agricultural Center (Parlier, CA, USA). This orchard was established in 1998 containing 51 verified hybrids. Each progeny genotype was represented by two trees in the orchard; the leaves were collected from parents and 51 Pop-DG populations, any one of the orchard tree and frozen in liquid nitrogen and stored at -80°C until used. Peach DNA was extracted from new leaves by the CTAB method (Doyle and Doyle 1987), modified by the addition of 1 % 2-mercaptoethanol to the isolation buffer. A preliminary quantification of DNA was determined after electrophoresis on 1 % agarose gels by comparison against a k HindIII standard of known concentration (Peace et al. 2005) (note: the GoldenGate Ò Genotyping assay requires using at least 15 ll of 50 ng/ll DNA per sample). The samples were sent to the DNA Technologies Core at the UC Davis Genome Center (http://dnatech.genomecenter. ucdavis.edu/dna_quant.html), and a second DNA quantification was done on a Molecular Devices Analyst plate reader with Pico Green (Invitrogen Molecular Probes). Phenotypic evaluation was carried out similar to our previous work (Ogundiwin et al. 2007). For genotyping of Pop-DG, five separate DNA samples of ''Dr. Davis'' and two ''Georgia Belle'' samples were used.

Discovery and selection of SNPs
The advantage of association mapping is that it can map quantitative traits with high resolution in a way that is statistically very powerful. Association mapping, however, also requires extensive knowledge of SNPs within the genome of the organism of interest, and is therefore difficult to perform in species that have not been well studied or do not have well-annotated genomes. In our present study, a set of 6,654 high-quality SNPs was developed under the framework of the ''Integration of genomic tools for next generation peach and almond cultivar development'' USDA-NRI project at UC Davis. SNPs were obtained by generation of whole genome sequence for ''Dr. Davis'', ''F8, 1-42'' and ''Georgia Belle'' using both Roche 454 and Ilumina-Solexa technologies (Ahmad et al. 2011). Assembly and alignment were done according to our previous published work (Ahmad et al. 2011). The SNPs were selected to be evenly distributed across the genome from a larger SNP set of 6,654 peach SNPs (approximately 1 SNP/ 40,000 nucleotides). SNPs for Pop-DG, a large number of these SNPs were found to be heterozygous for both parents, providing 1:2:1 segregation ratios. SNPs were also evaluated with a ''design'' score generated by Illumina, which attempts to predict the probability of success when used in a GoldenGate Ò assay. SNPs with scores = 1.1 were used for mapping. SNPs were named starting with ''UCD'' for the University of California, Davis, followed by ''SNP'' for the marker type and the number of their order in the OPA (UCD_SNP_XX). Once the final SNPs were selected, their reproducibilities and heritabilities were examined for Pop-DG population. Linkage group ( Oligonucleotide pool (OPA) development Oligonucleotides were designed, synthesized, and assembled into this OPA by Illumina Inc. The peach ''GS0012410-OPA.opa'' consisted of 1,536 SNPs and was used to genotype the Pop-DF and Pop-DG populations (Martínez-García et al. 2012). This is a unique set of SNP oligonucleotides developed for our previous project, but could be applied to any number of peach mapping/evaluation projects (Ahmad et al. 2011). For our present study, 1,536 SNPs genotyped from ''Dr. Davis'' and ''Georgia Belle'' were used. 960 SNPs with Minor Allele Frequency (MAF C10) were used for further analysis.

SNP genotyping
The SNPs were evaluated in an Illumina GoldenGate Ò Genotyping assay (Fan et al. 2003) with the iScan readout at the UC Davis Genome Center. The results were analyzed with Genome Studio TM Genotyping Module v1.0 from Illumina (Fig. 1). The SNPs were manually edited and removed from the analysis if clustering (segregation) errors in the parents/progeny distributions were observed, as described in the Infinium Ò Genotyping Data Analysis from Illumina Ò Technical Note (http://www.illumina.com), except that the Gen-Call score threshold level for inclusion of data into clusters was changed from 0.15 to 0.25 (Higher Gen Call scores indicate that the data points are more reliable and can be included in a specific cluster with more confidence). The deviations from expected allelic Mendelian inheritance ratios at all loci for each progeny were measured according to our previous paper (Martínez-García et al. 2012).

Association analysis
The population structure was inferred using the software program STRUCTURE 2.2 (Pritchard et al. 2000) on 960 SNPs (MAF C10 %) as described (Hao et al. 2012). In brief, five independent iterations of running were performed with hypothetic number of subpopulation (k) ranging from 1 to 10, based on burn-in period and the number of Markov chain Monte Carlo (MCMC) replications after burn-in, which were all assigned at 100,000. The correct estimation of k was provided by joining the log probability of data [Ln P(D)] from the STRUCTURE output and an ad hoc statistic Dk (Evanno et al. 2005), which was based on the rate of change in the log probability of data between successive k values. Based on k = 2 all 51 progenies were assigned subpopulation. The population structure matrix (Q) was generated for further analysis. To account for the population structure, four statistical models were evaluated: the general linear model (GLM) model without considering Q and K, the Q model with considering Q (not shown), the K model with considering K (not shown) and mixed linear model (MLM) model considering Q and K (Yu and Buckler 2006). Results of only two models, GLM model without considering Q and K and MLM model considering Q and K, were shown. Genome-wide association analyses based on these models were conducted with the software TASSEL 2.1. (Buckler et al. 2009) Markers were defined as being significantly associated with traits on the basis of their significant association threshold (-Log P C 2.00, P \ 0.0001, P \ 0.001, P \ 0.01 and P \ 0.1). (Buckler et al. 2009;Ge et al. 2003).

Results and discussion
Genome-wide association analysis (GWAS) GWAS have been more successful in crop plants than in humans (Brachi et al. 2011). GWAS overcomes several limitations of traditional gene mapping by (1) providing higher resolution, often to the gene level, and (2) using samples from previously well-studied populations in which commonly occurring genetic variations can be associated with phenotypic variation (Brachi et al. 2011). Using two models namely GLM and MLM, associations between CI traits and 960 SNPs genotypes were evaluated in our studied Pop-DG population. 22 SNPs were identified to have significant marker-trait associations at different threshold levels (-Log P C 2.00, P \ 0.0001, P \ 0.001, P \ 0.01 and P \ 0.1) for CI traits. The two models were compared based on different levels of significance across years at one or two time points. GLM analysis uses a R-restricted model that is a good test of marker effects, even when the data are unbalanced (Bradbury et al. 2007;Buckler et al. 2009). Similar to GLM, MLM performs an association test for each combination of traits and markers. 960 SNPs (MAF C10 %) were tested for association with CI traits like M, FB, and FBL. Results were observed and reported at four different significance levels with both models. Only SNPs with minimum significance level (P \ 0.1) at least 3 years were considered without considering the time points.

Mealiness (M)
Among 22 (Table 2). Higher proportion of significance (44 %) was found in MLM model (P \ 0.1) followed by GLM model (P \ 0.01). Trait FB did not show higher significance levels (P \ 0.0001) of association for all 12 SNPs (Fig. 1). More than one significant SNPs associated with QTL was identified for flesh browning in LG5 (qBrL5) in previous study (Martínez-García et al. 2012). CG1 was at same position of significant marker identified (UCD_SNP_1422) for flesh browning and CG38 was 4.5 cM away from qBrL5 (Dhanapal et al. 2012). In our present association analysis one SNP, UCD_SNP_1422, was found along with 11 others SNPs that showed association for browning.
Marker-trait association for CI Even though we used two different models in our present study, results from GLM model without considering Q and K and MLM model considering Q and K analysis did not show major difference in P values for all the traits under study. There was no additional markers pop out in either analysis when compared with one another but P values obtained from MLM analysis were lower than GLM (Tables 1, 2, 3). Two other models of GLM namely Q model with considering Q, the K model with considering K (results not shown) also resulted with same number of markers for association of various traits under study. Studies on soybean using various statistical models evaluated: the GLM model without considering Q and K; the Q model with considering Q; the K model with considering K; the MLM model with considering Q and K. MLM model performed a little better than the K model. So, the GWAS for soybean yield and yield components with the MLM model (Q ? K) was used to correct for population structure and genetic relatedness, using 1,142 SNPs and 209 haplotypes (Hao et al. 2012).
In plants, QTLs were originally mapped in biparental crosses, but they were restricted in allelic diversity and in having limited genomic resolution. The advent of highdensity SNP typing allowed whole-genome scans to identify often small haplotype blocks that are significantly correlated with quantitative trait variation. These approaches on plant studies have been successful in identifying loci that explain large portions of phenotypic variation, which is now available for most of the fruit tree species. The genetic control of CI related traits in peach has been studied and it has been demonstrated that mealiness, browning, and bleeding are probably controlled by major genes (Ogundiwin et al. 2007;Peace et al. 2006). Moreover, one major QTL has been detected for each of these symptoms of CI in LG4 and LG5, using a linkage map constructed from two segregating populations Pop-DG ('Dr. Davis' 9 'Georgia Belle') and Pop-G ('Georgia Belle' selfed) (Ogundiwin et al. 2007;Peace et al. 2006).
In our previous study a major QTL for mealiness and bleeding was found at the F-M locus at the bottom end of LG4. Other minor QTLs for mealiness were also found on LG4 and LG6. Besides, an Expressed Sequence Tags (ESTs) database has been developed specifically to study CI (Ogundiwin et al. 2008). Microarray analysis involving these ESTs has identified several cold-regulated peach genes some of which have been mapped close to CI QTLs on Pop-DG (Ogundiwin et al. 2009). A gene encoding a cell wall modifying enzyme, endopolygalacturonase (endoPG) co-localized with the major QTL affecting mealiness (Callahan et al. 2004;Peace et al. 2005). Another gene in the anthocyanin biosynthesis pathway, leucoanthocyanidin dioxygenase (PpLDOX), mapped to the same genomic region where the major QTL controlling browning was identified (Ogundiwin et al. 2009). Even though several previous studies were on CI traits using molecular markers SSRs, Restriction fragment length polymorphism (RFLP), and Single sequence conformational polymorphism (SSCP), this is the first study using SNP-based markers for CI traits.

Genetic-based association studies in Prunus
Association mapping (AM) is one of the effective approach for connecting phenotype and genotype in plants. It enhances previous QTL information for MAS in rice (Agrama et al. 2007), wheat (Tommasini et al. 2007), and maize (Yu and Buckler 2006). Development of genetic markers is crucial for association studies. A whole genome or selected segments of a genome of crop, plants, and trees are sequenced to identify differences across the genome. Subsequently, identified polymorphisms are genotyped across a larger and more diverse yet unrelated population. Prunus fruit development, growth, ripening, and senescence all include major biochemical and sensory changes in texture, color, and flavor. The genetic dissection of these complex processes has important applications in crop improvement (Cao et al. 2012). Fan et al. (2010) constructed a linkage map on peach F2 population of 378 genotypes developed from two genotypes with contrasting chilling requirements for QTL mapping. The study detected five QTLs of chilling requirement. Among these QTLs, qCR1a and qCR7 showed very prominent effects and were declared to be the major QTLs. The two QTLs not only facilitate marker-assisted breeding for low chill cultivars, but also paved the way for future fine mapping and mapbased cloning of genes controlling the chilling requirement. Another most highly significant association found was CPPCT005 on chromosome 4, which explained 25.1 % of the phenotypic variation (Fan et al. 2010).
Another study on apricot (Prunus armeniaca) for chilling injury requirement (CR) based on a two-way pseudo test cross mapping strategy, in which two highdensity apricot maps were constructed using a total of 43 SSRs and 994 amplified fragment length polymorphism (AFLP) markers that span an average of 502.6 cM with an average marker interval of 0.81 cM. Twelve putative CR QTLs were detected using composite interval mapping, a simultaneous multiple regression fit and an additive by-additive epistatic interaction model without dominance. An average of 62.3 ± 6.3 % of the total phenotypic variance was explained. QTLs corresponding to map positions of differentially expressed transcripts and candidate genes controlling a majority of the QTLs were shown to be stable between both Prunus species, as well as similar trends in their QTL effects, with the allele for increasing the trait value mostly originating from the high chill parents. The comparative QTL mapping strategy presented reveals the transferability of genetic information between two Prunus species, the characterization of stable QTLs, the utility of the maps to consolidate each other and to further validate previously identified CR QTL as a major controlling factor driving floral bud break.
Identifying the genetic variants that underlie complex traits was very essential to plant genetics. Two main approaches are available for mapping the relevant genes and identifying the variants associated with this complex trait: linkage mapping in families and population-based genetic association studies (Agrama et al. 2007). In theory, genetic association mapping is more powerful than linkage studies in identifying variants with weak effects that may contribute risks to common complex traits (Risch and Merikangas 1996). These trait-associated SNPs from our study were compared with other reported SNPs and SSR markers from QTL mapping analyses of our previous study on selected populations. Among associations, some were in regions where QTLs associated with the given trait had previously been identified. However, some associations were not consistent with the results of other published linkage maps because of different genetic materials and different marker systems used (Cao et al. 2012).

Conclusion
In our present study, association genetics was performed for CI susceptibility traits using GWAS based on GLM model without considering Q and K and MLM model considering Q and K analysis. A total of 22 SNPs markertrait associations were identified across 4 year in one or two time points by analyzing both models. However, further studies are necessary to confirm the association results in specific biparental crossing populations and consequently avoid identifying spurious associations. All three SNP markers namely UCD_SNP_1084, UCD_SNP_1441 and UCD_SNP_1172 will be investigated based on their position in genome, and it would be interesting to see if these SNPs are located in regulatory region or coding region in some candidate gene of interest. The presence of considerable levels of synteny among many rosaceous species suggests that even for those plant species with fewer genomic resources, candidate-gene association mapping coupled with QTL mapping studies and comparative mapping would be feasible and highly valuable. Knowledge acquired in one species can then be extended to others. Establishing a model species and refining its genome sequence and identifying regions associated with important traits could even be used for other related species of Rosaceae by conducting comparative genomics. Our future studies include localization of genes based on SNP marker-trait association for CI traits.