Association mapping for yield traits in Paeonia rockii based on SSR markers within transcription factors of comparative transcriptome

Background Allelic variation underlying the quantitative traits in plants is caused by the extremely complex regulation process. Tree peony originated in China is a peculiar ornamental, medicinal and oil woody plant. Paeonia rockii, one of tree peony species, is a precious emerging woody oil crop. However, in this valuable plant, the study of functional loci associated with yield traits has rarely been identified. Therefore, to explore the genetic architecture of 24 yield quantitative traits, the association mapping was first reported in 420 unrelated cultivated P. rockii individuals based on the next-generation sequencing (NGS) and single-molecule long-read sequencing (SMLRS). Results The developed 58 pairs of polymorphic expressed sequence tag-simple sequence repeat (EST-SSR) markers from 959 candidate transcription factors (TFs) associated with yield were used for genotyping the 420 P. rockii accessions. We observed a high level of genetic diversity (polymorphic information content, PIC = 0.514) and low linkage disequilibrium (LD) between EST-SSRs. Moreover, four subpopulations in the association population were revealed by STRUCTURE analyses. Further, single-marker association analysis identified 141 significant associations, involving 17 quantitative traits and 41 EST-SSRs. These loci were mainly from AP2, TCP, MYB, HSF, bHLH, GATA, and B3 gene families and showed a small proportion of the phenotypic variance (3.79 to 37.45%). Conclusions Our results summarize a valuable collection of functional loci associated with yield traits in P. rockii, and provide a precious resource that reveals allelic variation underlying quantitative traits in Paeonia and other woody oil crops.


Background
Woody oil crops are economically essential crops globally because of their high oil content in their seeds and/ or fruits, strong resistance, and stable yield. Now they have become the critical source for human edible oil, lubricants, biodiesel, cosmetics, paint and other industries [1,2]. P. rockii (S. G. Haw et L. A. Lauener) T. Hong et J. J. Li belongs to the Paeonia section Moutan DC., Paeoniaceae, which is originated from northwest China and has become one of the most representative species of tree peonies. The varieties originated from the species have been cultivated widely in China as the secondlargest cultivar group of Chinese tree peonies, including about 300 cultivars [3,4]. In addition to ornamental and medicinal cultivation, P. rockii also presents an advantage agronomic trait for its seed yield and unsaturated fatty acid contents in seeds in recent years [5]. Therefore, P. rockii has been rapidly extended as an emerging woody oil plant. However, how to improve the breeding efficiency and seed yield by using the molecular breeding methods has become a top priority. The breeding system of tree peonies is predominantly outcrossing and the cultivated cultivars are of hybrid origin [3]. Conventional breeding methods for tree peony germplasm mainly include cross and selective breeding, which takes at least 10 years to develop a stable new cultivar. Therefore, the long juvenile phase and complex genetic background make it very difficult to improve the yield using traditional breeding and reverse genetics methods. The better strategies of genetic architecture in the yield traits of P. rockii will be required to shorten the breeding cycle and improve the production value effectively.
Currently, association mapping has been widely used in the identification and genetic testing of important quantitative trait locus (QTL) in various species [6][7][8][9][10][11][12][13][14][15][16][17]. And some progress has been made in studies on crop yield. In popular Upland cotton (Gossypium hirsutum L.), 172 cultivars in China and 331 polymorphic SSRs were used for association mapping of yield-related traits. Totally, 93 significantly associated loci for seven yield traits were identified across more than one environment [18]. To reveal the genetic variations of yield and yield components traits in upland cotton, 403 accessions and 560 genome-wide SSRs were used for the association mapping based on a mixed linear model. A total of 43 marker loci were detected according to the best linear unbiased prediction and in at least three of the six environments (− lgP > 1.30, P < 0.05) [19]. In Hordeum vulgare, a total of 379 cultivars were used for genome-wide association mapping to identify the alleles controlling yield-related traits, and 13 putative genes regulating grain traits in European cultivated barley were obtained [20]. In wheat, association mapping was performed with a mixed linear model to identify the KASP marker for yield-related traits. The results showed that Hap-5A-1/2 of TaSnRK2.9-5A was significantly associated with high thousand kernel weight, while Hap-5A-4 with high grains per spike [21]. These provide a useful reference to identify the molecular markers which are closely associated with yield traits in woody oil crops including tree peonies.
Among dozens of markers, although SNP markers are usually used in crop association mapping, SNP sites are generally biallelic polymorphisms, which are significantly lower than SSR markers, and the development and detection costs are higher. Therefore, SSR markers have become increasingly popular in association analysis because of their codominant inheritance, extensive genome coverage, chromosome-specific location, relative abundance and high throughput genotyping [22]. In particular, the SSR loci in the coding regions, as the expressed TFs, are species conserved and can directly influence the gene transcription or translation [23][24][25]. Currently, EST-SSRs have used for the evaluation of genetic diversity and population structure, variety identification and map construction of different species [26][27][28][29][30][31][32][33]. And EST-SSRs have been considered to be the ideal marker type for the detection of woody plants with high heterozygosity. In tree peonies, reports on the development and application of SSRs are also increasing extensively [34][35][36][37][38][39][40]. However, the coverage of molecular markers, especially EST-SSRs, is still limited, which limits the extensive application of functional markers in genetic diversity, population structure, LD and association analysis.
TFs can modulate gene expression and play a essential critical role in the regulation of yield. For example, IPA1 (Ideal Plant Architecture 1) can promote both yield and disease resistance by sustaining a balance between growth and immunity in rice; The CCT domaincontaining gene family has significant impacts on heading date, regional adaptation and grain yield [41,42]. In wheat, TaNAC2-5A overexpressing transgenic wheat lines show higher grain yield and higher nitrogen accumulation in aerial parts; A wheat CCAAT Box-binding TF can increase the wheat yield with less fertilizer input [43,44]. In Maize, SBP-box TFs unbranched2 and un-branched3 can affect yield traits by regulating the rate of lateral primordia initiation; Nuclear factor Y (NF-Y) B subunits confer drought tolerance and lead to improved corn yields on water-limited acres [45,46]. In Brassica napus, overexpression of the brassinosteroid biosynthetic gene DWF4 simultaneously increases seed yield and stress tolerance [47]. In tree peonies, it is the first attempt to combine the SSR loci within TFs associated with yield to conduct the association mapping in P. rockii.
Taken together, the studies on the allelic variation associated with yield traits in P. rockii by combining association mapping and markers within TFs are rarely reported. Here, we sampled the cultivated P. rockii population for association mapping. After that, the genetic diversity and population structure were analyzed based on the developed EST-SSRs by comparative transcriptome. Moreover, LD and single-marker association mapping were conducted to explore the allelic effects on the natural variation of complex yield traits in the constructed population. Our results will lay a foundation to identify the linkage loci of yield traits in P. rockii and will be of great significance for the genetic improvement of yield traits of woody oil crops.

Distribution and statistics of phenotypic traits in the association population
In this study, multiple methods were used to conduct descriptive statistics on phenotypic traits and determine the traits affecting the yield of P. rockii. The one-way ANOVA analysis by R package showed a wide range of phenotypic variation among all the 24 quantitative traits measured. The variation coefficient ranged from 12.03 to 106.63% (mean 38.60%). The average variation coefficient of branch and leaf, flower and fruit traits was 22.57, 37.21 and 49.75%, respectively (Additional file 1: Table S1-1). Correlation analysis between different traits showed 202 significant correlations (P < 0.05), among which 182 showed a highly significant correlation (P < 0.01) (Additional file 1: Table S1-2). The principal component analysis showed that the 24 quantitative traits were divided into six principal components, and the cumulative variance reached up to 73.93% (Additional file 1: Table S1-3). Then, systematic cluster analysis was carried out based on the six principal components. The results showed that the association population was divided into eight categories and 24 quantitative traits were divided into five categories where the square Euclidean distance was respectively 12 and 18 (Additional file 1: Table S1-4, Additional file 2: Figure S1).
To sum up, when individual seed fresh weight was considered as the yield index, we found that each trait was related to yield. Simultaneously, individual fruit number, individual fruit fresh weight, and individual seed number could be considered as the primary traits affecting the yield of P. rockii. A total of 15 traits including flower diameter, petal length, petal width, plant height, crown breadth, compound leaf length, compound leaf width, maximum basal branch angle, single fruit length, single fruit width, single fruit pericarp thickness, multiple fruit fresh weight, multiple fruit seed number, multiple fruit seed fresh weight, and multiple fruit pod fresh weight could be considered as the secondary traits. The remaining traits could be deemed to beviewed as the auxiliary reference factors.

Genetic diversity
All the polymorphic 58 SSRs within TFs were used to evaluate the genetic diversity across 420 genotypes of cultivated P. rockii. Then, a total of 483 alleles were detected. The alleles per locus (N A ) ranged from 3 to 16, with an average of 8. The effective number of alleles (N E ) ranged from 1.02 to 7.32, with an average of 2.83.
The PIC values ranged from 0.021 to 0.849 (mean 0.514) and Shannon's information index (I) ranged from 0.07 to 2.18 (mean 1.13). The mean values of observed heterozygosity (Ho) and expected heterozygosity (H E ) were 0.529 and 0.561, respectively. Since Ho was higher than H E at 26 EST-SSR loci, the Wright's inbreeding coefficient (F IS ) ranged from − 0.828 to 0.867 with a mean value of 0.048 (Table 1). In conclusion, all the EST-SSRs showed a high level of polymorphism.
Population structure STRUCTURE analyses were used to determine the population structure. The result indicated that the Ln P(D) reached a mode at K = 4 before decreasing, and the highest delta K was detected when K = 4 ( Fig. 1a, b). The species genetic structure is dicussed with the results up to K = 6. Bayesian methods implemented in the STRUCTURE revealed extensive admixed ancestry for each sampled of P. rockii. Two major well-separated genetic clusters (I-III, IV) were identified at K = 2-6, among which four major well-separated genetic clusters (I, II, III, and IV) were observed at K = 4-6 ( Fig. 1c). Therefore, the 420 accessions were divided into four subpopulations.
Association mapping of yield quantitative-related traits in P. rockii Based on the genotype data, the Q matrix, the K matrix, and yield quantitative traits data, the MLM model was used to analyze the marker-trait associations. For floral traits associated with yield, we performed 232 (58 EST-SSRs × four traits) marker-trait association tests. Out of these, ten associations (4.31%) including 7 EST-SSR loci were significant at the P < 0.01 level. However, correction for multiple testing using the FDR method reduced the number to 8 (3.45%) at a significance threshold of Q < 0.05, involving three traits with 7 EST-SSRs. These loci explained 4.93 to 25.32% of the phenotypic variation, with an average of 11.49%. The number of associations varied across the traits was 2 (petal length) to 3 (flower diameter and petal number). One EST-SSR (PS31) was associated with at least one trait. For 3 out of the eight associations, the gene effect showed overdominance (|d/ a| > 1.25). The remaining 5 associations showed additive (|d/a| < 0.50, 2) or partial to full dominance (0.50 < |d/a| < 1.25, 3) (Additional file 3: Table S2-1).
For branch and leaf traits associated with yield, we performed 464 (58 EST-SSRs × eight traits) marker-trait association tests. Among them, 102 associations (21.98%) including 39 EST-SSR loci were significant at the threshold of P < 0.01. However, correction for multiple testing using the FDR method reduced the number to 87 (18.75%) at a significance threshold of Q < 0.05, involving seven traits with 36 EST-SSRs. These loci explained 3.79 to 37.45% (mean 11.06%) of the phenotypic variation. The number of significant associations varied across the traits was 2 (crown breadth and maximum For fruit traits associated with yield, we performed 696 (58 EST-SSRs × twelve traits) marker-trait association tests. Out of these, 64 associations (9.20%) including 30 EST-SSRs were significant at the threshold of P < 0.01. However, correction for multiple testing using the FDR method reduced the number to 46 (6.61%) at a significance threshold of Q < 0.05, involving seven traits with 26 EST-SSRs. These loci explained 4.45 to 30.29% (mean 11.47%) of the phenotypic variation. The number of significant associations varied across the traits was 3 (single fruit pericarp thickness, multiple fruit seed number) to 12 (single fruit number), as shown in Additional file 3. Eleven EST-SSRs showed significant associations with at least one trait. For example, PS43 was associated with single fruit number and multiple fruit pod fresh weigh. The modes of gene effect for 21 associations were overdominance, 15 were additive, and the remaining 10 were partial to full dominance (Additional file 3: Table S2-3).  Table S2). The gene families associated with all the three types of traits were AP2, GATA, HSF, and NAC, among which TFs from AP2 (29) gene family showed the highest association frequency. Moreover, all three types of traits exhibited the associations between identical traits and EST-SSRs from multiple gene families. For example, multiple fruit fresh weight was significantly associated with PS2 from the MADS-box gene family and PS30 from the AP2 gene family. Additionally, pleiotropism was generally shown in all three types of traits. For example, PS53 from the TCP gene family was significantly associated with multiple fruit seed number and multiple fruit seed fresh weight. PS12 from the MYB gene family was significantly associated with multiple fruit fresh weight, multiple fruit seed fresh weight, and multiple fruit pod fresh weight (Additional file 3: Table S2).

SSR loci within TFs associated with fruit traits
As fruit traits were significantly associated with yield, we further dissected the results of association mapping based on the EST-SSRs associated with fruit traits. Thus, five traits with small contribution rate (single fruit number, petal number, parietal lobule area, parietal lobule length and parietal lobule width) were excluded. Further, the results showed that 49 associations were significant at the threshold of P < 0.01 and Q < 0.05, which explained 4.45 to 37.45% (mean 11.09%) of the phenotypic variance. The number of significant associations varied across EST-SSRs was 1 to 7, and 11 out of the 19 EST-SSRs exhibited significant associations with at least one trait. The modes of gene effect for 21 associations were overdominance, 15 were additive, and the remaining 13 were partial to full dominance (Table 2). PS66 (HSF), PS2 (MADS-box), PS53, PS57 (TCP), PS7 (RING), PS145 (WD), PS122, PS131 (WRKY), PS94 (YABBY), PS85 (NAC) and PS43 (bHLH) were the EST-SSRs only associated with fruit traits, so these loci could be considered as the important references for MAS breeding in P. rockii.

Discussion
Genetic diversity of the population A vital prerequisite for dissecting population evolution and association analysis is to make clear the genetic diversity of the mapping population and the polymorphism of markers used [48]. Here, we evaluated the genetic diversity of the association population with 58 EST-SSRs developed in P. rockii. The results indicated that the average number of allele was 8, which was higher than that of Ficus carica [27], Pistacia atlantica [28], Cucumis melo [29] and cultivated P. rockii [39], but lower than that of wild P. rockii [4]. Here, we speculate that the differences should be due to the fact that the association population in this research is mainly composed of cultivated germplasm resources and there is a high level of heterozygosity in P. rockii. Simultaneously, the number and type of EST-SSRs and the characteristics of the test population are different.
Additionally, the population F IS [49,50] ranged from − 0.828 to 0.867 with an average of 0.048, of which 26 EST-SSRs were negative, indicating that there was significant heterozygous redundancy in the mapping population. We speculate that this is related to the hybrid origin and self-incompatibility of cultivated germplasm in P. rockii [3], which is consistent with the reported results in Populus tomentosa [51] and P. rockii [39], etc. Moreover, the positive selection of mutation loci and heterosis in the process of evolution is also critical reasons for heterozygous redundancy of P. rockii. Therefore, our results prove that highly heterozygous redundancy is an important biological characteristic of P. rockii, and further support the view that the cultivated resources in P. rockii are complex hybrid origin [3].
Comparing to other outcrossing woody plants, we detected a high level of genetic diversity in P. rockii (PIC = 0.514) [52,53], which was higher than that of P. tomentosa [51] and F. carica [27], but lower than that in P. atlantica [28] and Prunus avium [54]. This is different from the reported results of genic-SSRs with low polymorphism level, which may be due to the hybrid origin of cultivated P. rockii or the fact that we sample the accessions with significant genotype differences in this study.

Population structure
The analysis of crop population structure can not only reflect the gene exchange and affinity among individuals, but also be the premise of association mapping, which is conducive to improving the mapping efficiency and avoiding the emergence of false positives [55]. In this research, the association population was classified into four subpopulations based on STRUCTURE analyses, which is different from the previous report of three subpopulations in P. rockii [39]. We speculate that it is due to the differences in the EST-SSRs and population size. Moreover, the genetic information of the cultivated germplasm in P. rockii is mainly divided into four subgroups, which is consistent with the source of cultivated germplasm in this study, as shown in Fig. 3. We speculate that the four subpopulations may correspond to four gene pools and similarly reflect their geographic origins. It is also an essential reason for the rich genetic diversity.

The LD level of woody plants
The analysis of LD level in population is not only the genetic underpinnings for association mapping, but also provides a reference for the selection of appropriate association analysis strategies. Generally speaking, the LD level is low in cross-pollination plants, such as P. tomentosa [6,7,56], Eucalyptus [57], Gossypium hirsutum [58] and P. rockii [39]. The tree peony is an outcrossing species, thus we observed a low level of LD between EST-SSRs in the association population. We speculate that a significant amount of human intervention such as cross, controlled pollination, and germplasm exchange are also essential reasons of low LD levels in P. rockii. Moreover, the LD level of genic-SSR loci is not representative of the whole genome or the gene interval region level, which may also lead to the low LD levels [59]. Additionally, the lack of genomic information and the unknown genetic distance of EST-SSR loci in P. rockii limit the analysis for the decay of LD between markers, which remains to be studied further.
Association mapping QTL for yield quantitative traits SSR markers based on candidate genes have higher genetic effects in regulating the expression and function of genes related to quantitative traits [60]. In this study, 41 EST-SSRs and 141 associations related to yield traits were identified in the association population of P. rockii. Among them, SSR loci from MADS-box, AP2, MYB, and other gene families showed high genetic effects. Simultaneously, there were not only significant associations between one trait and SSRs in multiple gene families, but also pleiotropism or co-localized associations, which was consistent with the reported results [6,21,39]. We speculate that these pleiotropic associations are useful for discovering the important genomic regions and valuable in trait improvement by using MAS. In this research, we detected 176 combinations associated with three different types of traits. Still, this number was reduced to 141 after correction for multiple testing, which further improved the accuracy of the association results. Moreover, studies have shown that differences in population size and structure can cause variations in the results of association mapping. A typical association population should be a combination of multiple independent and unrelated individuals from the same region [61]. Therefore, to reduce the false positive association and provide a high precision estimation of allelic variation, we still need to verify the association results in combination with the validation populations in different regions or molecular biology experiments [62,63].
Association mapping of complex quantitative traits in plant can detect many significantly associated markers, but explain a small portion of phenotypic variance [64]. In this study, the average interpretation rate of EST-SSRs on flower (10.81%), branch and leaf (10.40%), and fruit traits (6.53%) was low, which were consistent with the reported results [59,61]. Besides, the analysis of genetic regulation relationship of quantitative traits in plants is helpful in utilizing the significantly association combinations for breeding. For example, the additive effect loci in the associations can account for significant genetic variation by cumulative genetic effects, while superdominant effect loci indicate that heterozygotes may superior to homozygotes. In this study, single fruit length was significantly associated with 2 EST-SSRs (PS131, PS94) with additive effects, which could add up to explain 14.63% of phenotypic variation without considering the intermarker effect. PS66 was significantly associated with  the multiple fruit seed fresh weight and showed a superdominant effect (|d/a| = 4.877), indicating that individuals with PS66 heterozygous loci in the association population might generate heavier seeds. So the combination of superdominant genetic loci might show greater heterosis. All in all, the combinations of genotypes with the same positive or negative effects can be used for the early selection of target traits.

TFs from AP2 gene family associated with yield
TFs with an APETELA2 (AP2) domain play significant roles in plant growth, development, and stress responses.
Studies have also shown that TFs from the AP2 gene family are associated with crop yield improvement [65][66][67][68]. In this research, TFs from the AP2 gene family were not only significantly associated with the three types of traits, but also showed the highest association frequency. Therefore, we speculate that TFs from the AP2 gene family are important factors in regulating yield quantitative-related traits of P. rockii. Moreover, TFs with an AP2 domain in Arabidopsis thaliana and Oryza sativa can control seed weight and seed yield [65,66]. In this research, PS19 (PB.50898.1) from AP2 gene family was significantly associated with seed number (multiple fruit seed number, |d/ a| = 1.517) and seed weight (multiple fruit seed fresh weight, |d/a| = 1.444) in P. rockii, and the modes of gene effect were all overdominance (|d/a| > 1.25). Therefore, we speculate that PS19 is likely to be an essential marker for regulating the seed yield in P. rockii. In addition, TFs from the AP2 gene family are also significantly associated with other yield traits such as fruit size and fruit weight, which remains to be further verified in future studies.

MYB-like TFs associated with yield
The MYB family of proteins comprise a large family of plant transcription factors, participating in the regulation of plant growth and development [69][70][71][72][73][74]. TFs with an MYB domain play an essential role in the regulation of crop yield [75][76][77]. In this research, in addition to the AP2 gene family, TFs from the MYB gene family also showed more associations in the mapping population. Three EST-SSRs (PS10, PS12 and PS17) from MYB gene family were significantly associated with fruit (4 associations) and branch and leaf traits (11 associations). Of these, PS12 (PB.59960.1) was significantly associated with fruit weight and seed weight-related traits, as shown in Additional file 3: Table S2. We speculate that PS12 is a key SSR loci for regulating yield of P. rockii. Also, a study has shown that the novel MYB-like TF OsMPH1 can regulate plant height and improve grain yield in rice [76]. We also observed that PS10 (PB.32979.2) and PS17 (PB.53756.4) from the MYB gene family were mainly associated with branch and leaf traits. Among them, only PS17 was significantly associated with plant height in the mapping population. So we speculate that PS17 is an important SSR loci regulating plant height associated with yield in P. rockii.

Conclusions
For the association mapping in P. rockii, we first constructed an association population consisting of 420 individuals. Then, the developed 58 polymorphic SSR loci within TFs associated with yield were selected for association mapping. Moreover, the genetic diversity and population structure were evaluated with the polymorphic loci, which proved that the population showed a high level of polymorphism and four subpopulations. Further, the results of association analysis based on single-marker showed that the 17 yield quantitative traits were regulated by 41 EST-SSR loci from 16 gene families, and 141 significant association combinations were identified. All these results furnish practical information to explore the effective functional loci associated with yield traits in tree peonies, it is also of great significance for the selection of yield-related traits and the cultivation of high-yield cultivars for oil tree peony.

Plant materials
More than 200, 000 seedling resources of cultivated P. rockii were mainly identified and collected from Gansu Province in northwest China by Beijing Guose Peony Technologies Co., Ltd. All these materials represented diverse genetic resources related to yield quantitative traits and were cultivated with general agronomic practices in the open field of Beijing Guose Peony Garden in Beijing, China (40°45′N, 115°97′E) (Fig. 3). Then based on the principle of covering the existing phenotypic variation as much as possible, we sampled a set of 420 individuals from the collection. All the evaluated samples were approximately 15 years old and covered three major flower types and eight color schemes. In total, 24 quantitative traits with stable performance were observed, demonstrating an effective assemblage of phenotypic traits.

Phenotypic data
The 420 individuals from the association population were scored based on 24  We have conducted an RNA-seq experiment of flower buds of P. rockii 'Jingshunfen' and P. rockii 'Fenmiantaosai' based on NGS and SMLRS technologies [79]. The RNA-seq data have been submitted to the NCBI Sequence Read Archive (SRR9915032, https://www.ncbi. nlm.nih.gov/sra/?term=SRR9915032, and SRR10872586, https://www.ncbi.nlm.nih.gov/sra/?term=SRR10872586). Based on the sequencing data, the candidate 959 TFs from 21 gene families associated with yield based on previous reports were screened. Then a total of 166 EST-SSRs containing six nucleotide repeat types were identified, with an average of one SSR per 5.78 unigenes. Among them, 58 polymorphic EST-SSR markers have been identified and proven to be effective in Paeonia (Unpublished observations 1 ). These 58 EST-SSRs were selected to genotype the 420 accessions (Additional file 5: Table S4).
Each PCR reaction was performed in a total reaction mixture volume of 10 μL containing 5 μL 2 × Power Taq PCR MasterMix (BioTeke, Beijing, China), 0.5 μL (10 pmol) each primer, 3.0 μL dd H 2 O and 1 μL (50 ng) template DNA. The amplification program was as follows: 5 min at 95°C, 35 cycles of 30 s at 95°C, 30 s at the appropriate annealing temperature, 1 min at 72°C, and 10 min at 72°C, 4°C hold. The products were separated by capillary electrophoresis using an ABI3730XL capillary sequencer along with an internal size standard (Applied Biosystems, Carlsbad, CA, USA) after confirmation of PCR amplification by electrophoresis on a 1% agarose gel. The polymorphic EST-SSR loci were read with Gen-eMarker v1.80 software using LIZ 600 size standards (SoftGenetics, State College, Pennsylvania, USA). Subsequently, Micro-Checker v2.2.3 (http://www.microchecker.hull.ac.uk/) was applied to identify and to correct the genotyping errors [80].

Data analysis
Descriptive statistics were performed for 24 quantitative traits, including coefficient of variation (CV/ % = standard deviation / mean × 100%), one-way ANOVA (single factor completely randomized trial design), correlation analysis, principal component analysis (PCA), and cluster analysis. All the statistics were carried out by using IBM SPSS Statistics 20.0 software and R language package.
The developed 58 polymorphic EST-SSRs were used to analyze the genetic diversity of association population. The summary statistics of N A , N E , Ho, H E , PIC, I and F IS were calculated by GenAlEx v6.501 [81] and POPGENE v1.32 [82]. The developed 58 polymorphic EST-SSRs were also used to analyze and evaluate the population structure. The Bayesian method in the software package, STRUCTURE 2.3.4 (http://pritch.bsd.uchicago.edu/ structure.html) [83], was used to infer the number of subpopulations (K) through an admixture model. For each value of K (K = 1-19), ten independent runs were performed with a burn-in period of 100,000 followed by 200,000 Markov Chain Monte Carlo (MCMC) replications. Then the results were submitted to the Structure Harvester (http://taylor0.biology.ucla.edu/struct_harvest/ ) [84]. As a result, LnP(D) and ΔK were used to detect the optimum K value [85]. Then the CLUMPP v1.1.275 was used to analyze the results from replicate analyses for optimal alignments of replicate clusters [86]. The output from CLUMPP was displayed by the cluster visualization program DISTRUCT [87]. Structural analysis was used to determine the optimum population structure and correct false positives for association mapping.
LD was measured as the squared correlation of allele frequencies r 2 . The r 2 values between pairs of EST-SSRs (minor allele frequencies > 1%) were calculated with 10 5 permutations using the TASSEL v2.0.1 software (http:// www.maizegenetics.net/). The pairs of loci were considered to show a significant LD when P < 0.001. In the association mapping, the TASSEL v2.0.1 software package was used for marker-trait analysis with 10 4 permutations by using the mixed linear model (MLM) [88]. The Q matrix estimating the membership coefficients for each accession was derived from the STRUCTURE runs. The relative kinship matrix (K) was determined by using SPAGeDi software v1.2 [89]. False discovery rate (FDR) analyses were conducted using QVALUE in R [90].
The gene effect was assessed by using the ratio of dominance (d) to additive (a) effects, which was estimated from least-square means for each genotypic class. The values of 0.50 < |d/a| < 1.25, |d/a| ≤ 0.5 and |d/a| > 1.25 were defined as partial or complete dominance, additive effects and under-or overdominance, respectively. The algorithms of dominance and additive effects were d = G Bb -0.5 (G BB + G bb ) and 2a = |G BB -G bb |, respectively (Gij: The mean value of the phenotype corresponding to the ijth genotype; BB, bb: Different homozygous genotypes; Bb: Heterozygous genotype) [91].