Breeding programs require exhaustive phenotyping of germplasms, which is time-demanding and expensive. Genomic prediction based on next-generation sequencing techniques helps breeders harness the diversity of any collection to bypass phenotyping. Here, we examined the genomic prediction’s potential for seed yield and nine agronomic traits using 26171 single nucleotide polymorphism (SNP) markers in a set of 337 flax (Linum usitatissimum L.) germplasm, phenotyped in five environments. We evaluated 14 prediction models and several factors affecting predictive ability based on cross-validation schemes. Most models gave close predictive ability values across traits for the whole marker set. Models covering non-additive effects yielded better predictive ability for low heritable traits, though no single model worked best across all traits. Marker subsets based on linkage disequilibrium decay distance gave similar predictive abilities to the whole marker set, but for randomly selected markers, it reached a plateau above 3000 markers. Markers having significant association with traits improved predictive abilities compared to the whole marker set, when marker selection was made on the whole population instead of the training set indicating a clear overfitting. The correction for population structure did not increase predictive abilities compared to the whole collection. However, stratified sampling by picking representative genotypes from each cluster improved predictive abilities. The indirect predictive ability for a trait was proportionate to its correlation with other traits. These results will help breeders to select the best models, optimum marker set, and suitable genotype set to perform an indirect selection for quantitative traits in this diverse flax germplasm collection.
Article
Genomic prediction for agronomic traits in a diverse Flax (Linum usitatissimum L.) germplasm collection
https://doi.org/10.21203/rs.3.rs-3213823/v1
This work is licensed under a CC BY 4.0 License
Journal Publication
published 08 Feb, 2024
You are reading this latest preprint version
Genome-wide association study (GWAS)
genomic prediction (GP)
canola
sclerotinia stem rot (SSR)
single nucleotide polymorphism (SNP)
marker-trait-associations (MTAs)
Flax (Linum usitatissimum L.), a natural source of oil and fiber, is grown throughout the world since prehistoric times and has considerable economic importance 1. Flaxseed plays an important role in human nutrition 2,3 by providing oil rich in omega-3 fatty acid, dietary fibers, and anti-carcinogenic lignans. Flaxseed oil is also used for different industrial purposes such as ink making, varnishing, painting, and road carpeting due to its specific drying properties 4. Additionally, flaxseed meal has value as poultry and animal feeds 5,6 while flax fiber is used in making linen cloth and different bio-industrial products 7.
Among the flax-growing countries, Kazakhstan produces the most oilseed flax followed by the Russian Federation, Canada, China, and the USA, whereas France alone produces three-fourths of fiber flax in the world followed by Belgium, Belarus, Russian Federation, and China 8. In the United States, North Dakota (ND) has the greatest % of cultivated flax acres, which covers about 71% (215 million acres) and 80% (3.75 billion bushels) of U.S. flax acreage and production, respectively, and annually contributes about $46 million U.S. dollars to the national economy (Data averaged from 2017 to 2021) 9. Like many other cultivated crops, flax production in ND is being challenged by different biotic and abiotic stresses 10,11. To combat these challenges and meet farmers desires for high-yielding varieties with increased oil and protein content, North Dakota State University (NDSU) runs a moderate-size flax breeding program. The program utilizes classical breeding methods, especially modified bulk methods to develop varieties, which is expensive, laborious, and time-consuming. To speed up the breeding process, it is prime time to adopt cutting-edge breeding tools such as marker-assisted selection (MAS) and genomic selection (GS) in the program.
Initially, breeders utilized marker-trait associations revealed by linkage mapping and genome-wide association mapping for MAS in breeding programs to enhance efficiency and genetic gain 12. To date, MAS has been successfully used to improve monogenic or oligogenic traits in many major crops such as rice 13,14, wheat 15–17, maize 18–20 etc. However, the improvement of quantitative traits controlled by multiple QTLs with minor effects is challenging. A multi-marker MAS system can be used to improve quantitative traits, but it is very difficult to identify and account for all the allele effects 21,22. Breeders can overcome the limitations of MAS by using a genome-wide selection approach. Genomic selection (GS), also known as genomic prediction, by considering all marker effects regardless of significance holds a promise to accelerate the rate of genetic gain in the case of quantitative traits 23,24. The GS was first successfully used in animal breeding, where it was applied to dairy cattle 25. In recent days, the low genotypic cost compared to phenotypic cost has made the GS an attractive decision tool to select and evaluate accessions in diverse germplasm collections.
To date, many statistical models for GS have been developed. Initially, the RRBLUP model was widely used 23,26. Later, plethora of linear or parametric models 27–30 and non-linear or non-parametric models 31–33 have been evolved. Linear models capture only additive gene action effects, whereas non-linear models capture additive and non-additive (dominance, epistasis, and pleiotropy) interactions. Different models vary in their underlying assumptions and algorithms and many authors confirmed that no single model worked best across traits, rather a particular model outperformed a particular trait 34–37. That is why it is always recommended to test multiple models across traits to achieve maximum prediction accuracy. Along with models, various factors such as relatedness between the training set and validation set 38,39, correlation among studied traits 40,41, trait heritability, marker density, QTL size with effects 34,42−46 and genetic diversity of studied collection 36,47,48 also affects the prediction accuracy.
Plant breeders have successfully utilized GS to accelerate the varietal development process in rice 49–52, wheat 53–57, maize 58–61 and other crops. Despite the wide application of GS in many crops, its utilization in flax breeding has not flourished yet. Because very few reports are available regarding the utilization of GS in bi-parental 62,63 and diverse 64 flax populations, there is great potential for evaluating GS in diverse flax collections. No research has been conducted before on GS targeting complex agronomic traits for the genetic exploration of the NDSU flax germplasm collection. In the present research, we used a collection of 337 flax germplasms and 26,171 single nucleotide polymorphism (SNP) markers to explore the potential of genomic prediction implementation for different agronomic traits, including determining the best prediction models across traits and optimum marker numbers. We also examined the impact of taking into account marker-trait association, population structure, and trait correlations in the genomic prediction framework on prediction ability.
Plant materials
We collected 500 flax accessions and their wild relatives from North Central Regional Plant Introduction Station (NCRPIS), Ames, Iowa, USA. All genotypes were grown in the field as single rows. We discarded the heterogeneous rows and kept the homogeneous lines for parental stock. Finally, we made a core collection of 337 flax germplasm accessions, which comprises homogeneous lines from NCRPIS, NDSU-released varieties and advanced breeding lines, and varieties developed by different institutes in the USA and Canada (Supplementary Table S1). The advanced breeding lines (F7 generation) were obtained by crossing different parents in various combinations. The core collection is being maintained through selfing.
DNA extraction, sequencing, and SNP calling
Young leaves collected from 30-day old plants were used as the source of DNA. The collected leaf samples from each genotype were lyophilized and subsequently pulverized using stainless beads in a plate shaker. DNA was extracted from the ground leaf tissue using a Qiagen DNeasy Kit (Qiagen, CA, USA) according to the manufacturer’s protocol. A NanoDrop 2000/2000c Spectrophotometer (Thermofisher Scientific) was used to measure DNA concentrations and for the genotype by sequencing (GBS) approach. The GBS library was prepared using ApekI enzyme 65 and sequencing of the library was accomplished using an Illumina HiSeq 2500 sequencer at the University of Texas Southern Medical Center, Dallas, Texas, USA. Identification of SNPs was based on a 120-base kmer length and minimum kmer count of ten using the TASSEL 5 GBSv2 pipeline 66 and the flax reference genome 67 (available at: https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/224/295/GCA_000224295.2_ASM22429v2/). The reads were aligned to the reference genome using Bowtie 2 (version 2.3.0) 68 alignment tool, which identified 243,040 SNPs that passed all required steps of the TASSEL 5 GBSv2 pipeline. Though flax is a strictly self-pollinating crop and inbred lines were used, there was a possibility of heterozygous SNPs due to artefactual collapse of homologous sites during alignment. We removed the heterozygous SNPs and filtered the row SNP set using VCFtools 69 following the criteria: minor allele frequency (MAF) ≥ 0.05, missing values (max-missing) ≤ 25%, depth (minDP) ≥ 3, min-alleles = 2 and max-alleles = 2. This filtering process yielded 26,171 bi-allelic high-quality SNP markers.
Phenotyping
We planted 337 genotypes following augmented row-column design 70 with three standard checks (ND Hammond, Gold ND, and Omega), and the checks were diagonally placed to cover spatial heterogeneity (Supplementary Fig. S1). Each check was replicated 20 times per trial and experiments were conducted at Fargo, ND, USA (46.8772° N, 96.7898° W) in three consecutive years (2018, 2019, and 2020) and at Carrington, ND, USA (47.4497° N, 99.1262° W) in two consecutive years (2019, 2020). Hereafter, we referred to the location-year combinations as environments: E1 (Fargo, 2018), E2 (Fargo, 2019), E3 (Carrington, 2019), E4 (Fargo, 2020), and E5 (Carrington, 2020). Within environments, for each genotype, we used 8 (4m×2m) m2 four-row plots and 30 gm of seeds per plot for planting except E1. In the case of E1, we used a single row per genotype, 7 gm of seeds per row. Standard fertilization and cultural practices were used throughout the experiment. Data for nine agronomic traits was recorded in all environments and seed yield in four environments. For each genotype, data measured on different traits was based on the criteria and methods described by Nôžková et al. (2011) 71 with minor modifications. Days to flowering was determined as the number of days from planting to when approximately 50% of plants per plot start flowering. Plant height was measured as the length of the main stem at maturity from the hypocotyl ending point to the plant’s top. Technical length was measured as the length of the main stem at maturity from the end of the hypocotyl to the point where branching starts. The branch number per plant was counted as the lateral branches of the main stem inflorescence. In this case, only the primary lateral branches were considered. Boll number per plant was the capsule number carried by the main stem inflorescence. Thousand seed weight was the weight of exactly 1000 seeds. Seed area, seed width, and seed length were calculated as an average of 1000 seeds. We used a MARVIN seed analyzer (GTA Sensorik GmbH) to measure seed-related traits. We harvested each plot separately and measured the grain weight in grams as yield per plot.
Phenotypic data analysis
In this study, a two-stage analysis of phenotypic data was performed. In stage-I, the best linear unbiased estimates (BLUEs) and other statistics for all genotypes within each environment were determined using the following model:
$$y=X\tau +{e}_{R}+{e}_{C}$$
1
Where X is the design matrix, τ is the fixed effect of genotype, eR is the random effect of row and eC is the random effect of the column.
In stage II, we fitted the BLUEs and weights from stage-I analysis in equation-2 and estimated the best linear unbiased predictions (BLUPs) of genotypes across all environments.
$${y}_{ij}=\mu +{G}_{i}+{E}_{j}+{GE}_{ij}+{e}_{ij}$$
2
Where \({y}_{ij}\)is the observed phenotypic value of the ith genotype in the jth environment, µ is the overall mean, \({G}_{i}\)is the random effect of the ith genotype, \({E}_{j}\) is the fixed effect of the jth environment, \({GE}_{ij}\) is the G×E interaction term and \({e}_{ij}\) is the residual error. The analysis was done using the R-shiny app MrBean (https://beanteam.shinyapps.io/MrBean/).
We calculated the heritability of each trait using the following formula proposed by Cullis et al. (2006) 72:
$${H}_{Cullis}^{2}=1-\left(\frac{PEV}{md\text{*}{V}_{g}}\right)$$
3
Where the genotypic predicted error variance is \(PEV\), \({V}_{g}\) is the genotypic variance and \(md\) is mean values from the diagonal of the relationship matrix. The heritability calculation was done using the R package Sommer 73.
We calculated Pearson correlation among different traits within each environment and combined all environments using observed unadjusted phenotypic values. The correlation of phenotypic values of a trait observed in different environments were also calculated; for this purpose, we used the R package corrplot 74.
Structure and linkage disequilibrium (LD) analysis
An admixture model-based structure analysis of the whole germplasm set was conducted using STRUCTURE 75 software utilizing the whole SNP marker set (26,171). To strengthen the result, structure analysis was run at various combinations of burn-in lengths (5000 to 50000) and Monte Carlo Markov Chain (MCMC) lengths (5000 to 100000). Each combination was replicated 10 times per K (K1 to K10). Along with the DeltaK approach 76, four alternative statistics (MedMedK, MedMeaK, MaxMedK, and MaxMeaK) 77 were used to identify optimum cluster numbers – StructureSelector 78 was used for this purpose. We assembled 10 replicates of the Q-matrix for the best-fitted cluster number using CLUMPP 79. Principal component analysis (PCA) was run using a covariance standardized approach in TASSEL 80. To show the genetic divergence among identified clusters, pairwise Fst81 was calculated using Arlequin3.5 82 at 10000 permutations. Gaussian finite mixture model-based clustering of the collection was fitted via the EM algorithm in R package mclust 83 using phenotypic data. We also visualized the PCA and phenotype-based clustering output using the ggplot2 R package 84.
Chromosome-wise LD (r2 values) among SNPs was calculated using the 26,171 SNP markers in PopLDdecay software 85. The LD decay rate was defined as a half-decay distance, at which observed r2 between sites decays to less than half of the maximum r2 value 86. For this purpose, we wrote R scripts (available on personal communication).
Genomic prediction models’ comparison
In the case of genomic prediction, the linear model equation is unsolvable as the explanatory variables (marker number) exceed the observation number. Researchers solve this problem by using ridge regression or Bayesian computations or parametric method and penalized regression or semi-parametric method 87. Generally, all methods are fitted to the basic skeleton (Eq. 4) with modifications:
Where y is the phenotypic value (BLUP), µ is the fixed intercept, Z is the marker matrix, u is the marker effect vector and ɛ is a residual vector.
In this study, we calculated the prediction ability for different traits using 26,711 markers and BLUP values based on 14 different parametric and semi-parametric models such as GBLUP 88–90, EGBLUP 91, RR 23,26,92, LASSO 93, EN 94, BRR 95, BA 23, BB 96, BC 97, BL 29, RKHS 98, RF 99, SVM 100,101 and MKRKHS 33. Details of the models are available in previously published research articles 34,95. For this purpose, we used the R package BWGS pipeline 34. For each model, we calculated the predictive ability as Pearson correlation between genomic estimated breeding values (GEBVs) and phenotypic values (BLUPs) of the validation set (VS). For this, a fivefold cross-validation approach was used, i.e., we randomly selected 80% of the collection as a training set (TS) and the remaining 20% as a validation set (VS). The process was repeated 100 times and finally predictive ability was reported as average across 100 replicates. The model that gave the best predictive ability for a particular trait was declared as the best-fitted model for the corresponding trait. For subsequent analyses, we only used the best-fitted model identified in this stage for specific traits.
Optimal marker density determination
Predictive ability for each trait was calculated using different subsets of the whole marker set. The marker subset was based on linkage disequilibrium (LD) decay distance and random selection and the marker subset was thinned using chromosome-wise LD decay distance, which yielded 5362 markers. We also selected subsets of 20, 200, 1000, 3000, 7000, and 13000 markers based on random sampling to minimize or avoid any biases. Using a five-fold cross-validation approach, predictive ability in each subset was measured. The process was repeated 100 times and finally, predictive ability was reported as average across 100 replicates.
Marker selection based on marker-trait association
Significant markers for each trait were grouped based on genome-wide association mapping (GWAS) results. This was done five different ways. In the case of scenario-I, we conducted SNP-based GWAS for all traits within each environment (BLUEs) and combined all environments (BLUPs) using different single locus and multi-locus models such as general linear model (GLM) 102, mixed linear model (MLM) 103, compressed MLM (CMLM) 104, enriched compressed MLM (ECMLM) 105, settlement of MLM under progressively exclusive relationship (SUPER) 106, multiple loci MLM (MLMM) 107, fixed and random model circulating probability unification (FarmCPU) 108 and Bayesian information and linkage-disequilibrium iteratively nested keyway (BLINK) 109. The R package GAPIT (version 3) 110 was used to run GWAS and the best-fitted model for a particular trait was determined based on the mean of squared difference (MSD) values and QQ plots. Using the best-fitted model output, we identified significant SNPs associated with a particular trait based on a p-value threshold. The p-value threshold was calculated by dividing the type-I error rate (α) by the effective number of independent tests (Meff) at α = 0.05 111. In this study, the p-value threshold was 0.000103. Any marker having a p-value ≤ p-value threshold was declared as a significant marker. We grouped the significant markers for each trait considering each environment (BLUEs) and combined environment (BLUPs). In the case of scenario II, GWAS was conducted using only an MLM model based on combined environment BLUP values. The top 20 markers for each trait having the lowest p-value were grouped. Then the selected marker set from scenario-I & II were used to calculate predictive ability 100 times using a five-fold cross-validation approach and predictive ability was reported as average across 100 replicates. In the case of scenario III, we randomly divided the whole collection into TS (80% of the whole collection) and VS (20% of the whole collection) 15 times. Each time we did GWAS using TS only and grouped the markers as described in scenario II. Then, predictive ability was calculated and reported as the average across 15 replicates. In the case of scenario IV, GWAS was conducted as one-way ANOVA using the R function lm, and every marker was tested one at a time using phenotype (BLUP values) considering the whole collection. The markers having a p-value less than the predefined p-value (0.001, 0.01, and 0.05) were used for predictive ability calculation 100 times with five-fold cross-validation. In the case of scenario V, we followed scenario III and used one-way ANOVA instead of MLM model-based GWAS to select the marker set having a p-value less than the predefined p-value. The selected marker set was then used to calculate predictive ability 100 times using five-fold cross-validation and predictive ability was reported as average across 100 replicates.
Prediction ability considering population structure
To investigate confounding effects of population structure on predictive ability, we calculated predictive ability using a subset of the whole collection after discarding the genotypes belonging to the most divergent (clusters showing the highest pairwise Fst value to other clusters) genetic clusters. Divergent clusters were removed to increase the genetic similarities among the remaining individuals. We incorporated the Q-matrix of structure analysis output into the RR-BLUP model to nullify the effect of structure and calculated predictive ability considering the whole collection. Albeit a smaller sample size, cluster-wise predictive ability. In all cases, we used five-folds cross-validation 100 times. To cover the genetic variance, we also did stratified sampling and calculated predictive ability following two different methods (M-I and M-II). The whole collection was arranged into small groups of 25, 50, 75, 100, 125,150, 175, 200, 250, and 300 genotypes by randomly selecting genotypes from each cluster proportional to the size of the cluster. In the case of method-I (M-I), 100 times five-fold cross-validation within each subset was done to report the predictive ability as the average across all runs. In the case of method -II (M-II), each subset was used as TS and the remaining genotypes as VS. For this purpose, we made 20 replicates of each subset. The predictive ability was calculated as Pearson correlation between GEBVs and BLUP values of VS and was reported as an average across 20 times. There was overlapping of genotypes among replicates as genotypes were selected randomly from each cluster proportional to the size of the cluster each time.
Indirect predictive ability calculation
Indirect prediction ability for each trait, considering other traits separately, was calculated using the five-fold cross-validation schemes 100 times. For example, we estimated the GEBVs of the validation set for plant height. Then we calculated the indirect prediction ability for seed yield by plant height as Pearson correlation between GEBVs of the validation set based on plant height and BLUP values of that set considering seed yield. The same was repeated 100 times and prediction ability was reported as average across 100 replicates. The same procedure was followed for different trait combinations.
Phenotypic variability
The genotype collection showed continuous variation for all traits under all environments (Supplementary Fig. S2). Under the five environmental conditions, days to flowering ranged from 36 to 67, plant height ranged from 20.0 to 82.8 cm, technical length ranged from 6.5 to 60.7 cm, branch number ranged from 2.7 to 12.3 and boll number ranged from 5 to 50. All these traits had variable coefficient of variation (CV) values across environments. Among these, for boll number, a maximum CV value of 36.2% was observed in E3. However, seed-related traits had relatively lower CV values compared to other traits. Among these, thousand seed weight had a maximum CV value of 16.0% in E1. Seed yield was evaluated under four environments, which ranged from 87 to 630 g per plot with a CV of 35.7 to 41.4%. We calculated environment-specific heritability as well as combined heritability considering all environments. Boll number and seed yield had the lowest heritability (< 0.40) compared to other traits under all environments. Days to flowering, branch number, and seed width indicated both low to high (0.18 to 0.93) environment-specific heritability, whereas it was high (> 0.60) for plant height, technical length, thousand seed weight, seed area, and seed length. All the traits had high combined heritability (> 0.65) except the branch number. The details of phenotypic variability are presented in Supplementary Table S2.
Phenotypic correlation
We investigated the phenotypic association among different traits and seed yield (Supplementary Fig. S3) in all environments. Days to flowering produced both positive and negative weak correlations with other traits across environments. Plant height, technical length, branch number, and boll number had a positive significant correlation among them in all environments except E1. Plant height and technical length were always highly correlated (r > 0.77) in all environments. Seed-attributing traits such as thousand seed weight, seed area, width, and length were significantly positively correlated with each other, but they exhibited mostly weak correlations with the remaining traits and seed yield across environments. Seed yield had a positive significant association with plant height, technical length, branch number, and boll number, but a weak positive and negative association with days to flowering in all environments. Among all the traits, the boll number had the best positive correlation (r = 0.69 to 0.79), and the seed area had the best negative correlation (r = -0.32 to -0.18) to seed yield.
We also calculated the correlation among environments for each trait (Supplementary Fig. S4) and found strong associations (r > 0.5) among environments for plant height, technical length, thousand seed weight, seed area, and seed yield. The same scenario was observed for boll number, seed width, and seed length whereas E5 had moderate correlations (r < 0.5) with other environments. For branch number and days to flowering, environments show comparatively low correlations among them.
Population structure and linkage disequilibrium (LD) analysis
Structure analysis revealed 3 to 9 clusters based on the Delta K approach 76 and 4 to 5 clusters based on four alternative statistics (MedMedK, MedMeaK, MaxMedK, and MaxMeaK) 77 (Table 1). Based on the structure output, we separated the whole collection into five clusters (P1 to P5). Winter-type Hungarian (European), spring-type Asian (Indian and Pakistani), and winter-type Turkish genotypes were dominant in clusters P1, P3, and P4, respectively. Cluster P2 contained mixed-type genotypes of different origins while cluster P5 was dominated by spring type NDSU advanced breeding lines of American origin and fiber type (Supplementary Table S1, Supplementary Fig. S5).
| Structure run # | Burn-in lengths | MCMC lengths | Number of clusters (K) | Number of Reps | Number of clusters | ||||
|---|---|---|---|---|---|---|---|---|---|
| ∆Kα | Med MedKβ | Med MeaKβ | Max MedKβ | Max MeaKβ | |||||
| 1 | 5000 | 5000 | 10 | 10 | 9 | 4 | 4 | 5 | 5 |
| 2 | 10000 | 10000 | 10 | 10 | 9 | 4 | 4 | 5 | 5 |
| 3 | 10000 | 50000 | 10 | 10 | 3 | 4 | 4 | 5 | 5 |
| 4 | 20000 | 20000 | 10 | 10 | 9 | 4 | 4 | 5 | 5 |
| 5 | 20000 | 50000 | 10 | 10 | 6 | 4 | 4 | 4 | 4 |
| 6 | 50000 | 50000 | 10 | 10 | 7 | 4 | 4 | 5 | 5 |
| 7 | 50000 | 100000 | 10 | 10 | 6 | 4 | 4 | 5 | 5 |
We also performed principal component analysis (PCA) to show the genetic similarity among genotypes. The first two axes explained 18% of the total observed variation. The output of principal component analysis was in line with that of structure analysis (Fig. 1A). A pairwise Fst comparison indicated genetic divergence among clusters. All combinations showed significant pairwise Fst comparison at p < 0.01. We found Fst ≥ 0.20 for all combinations except combinations P1 and P2, P2, and P4. Cluster P3 showed maximum divergence (Fst ≥ 0.30) from all other clusters (Table 2).
| Cluster pairwise Fst | |||||
|---|---|---|---|---|---|
| P1 | P2 | P3 | P4 | P5 | |
| P1 | 0 | ||||
| P2 | 0.13** | 0 | |||
| P3 | 0.48** | 0.38** | 0 | ||
| P4 | 0.21** | 0.13** | 0.47** | 0 | |
| P5 | 0.24** | 0.20** | 0.50** | 0.30** | 0 |
Moreover, Gaussian finite mixture model-based clustering of the whole collection using phenotypic data yielded four clusters (C1 to C4) (Fig. 1B). Cluster C1, C3, and C4 contained genotypes of different types and origins, while cluster C2 contained the highest number of genotypes dominated by spring type NDSU advanced breeding lines (Supplementary Table S1). Phenotype-based clustering was not consistent with genotype-based clustering (Fig. 1C, Fig. 1D) i.e., genotypic clusters containing genotypes belong to different phenotypic clusters and vice-versa.
In the whole collection, LD decayed to its half maximum within < 21 kb. LD decay rate varied according to chromosome (Supplementary Fig. S6), which was lowest in chromosomes Lu1 and Lu3 (32 kb) but was highest in chromosomes Lu7 and Lu8 (15 kb).
The efficiency of different genomic prediction models
We ran 14 different genomic prediction models available in the BWGS package to estimate the prediction ability in 100 independent five-fold cross-validations using 26,171 markers for ten agronomic traits. We determined the efficiency of each model in terms of time requirement and predictive ability (Fig. 2 and Fig. 3). For all traits, GBLUP required less time (< 60 min.) except the trait thousand seed weight (567 min.) and RF was the most time demanding (1637 to 2414 min.) model. EGBLUP, a modification of the GBLUP model, which covers epistatic interactions, required more time (767 to 1428 min.) than GBLUP. Besides GBLUP, other less time-demanding models were EN (51 to 124 min.), LASSO (82 to 143 minutes), SVM (120 to 155 min.), and RR (233 to 374 min.). As time requirement by different models is heavily affected by computer configurations, we always used the same set-up for each run and the values presented here must be taken only for comparison.
The predictive abilities of different models vary according to traits (Fig. 3, Supplementary Fig. S7). For all traits, SVM yielded the poorest predictive ability (-0.03 to 0.09). However, SVM, LASSO and EN also resulted in low predictive ability for all traits. Model RR yielded the highest predictive ability for plant height (0.48), technical length (0.56), boll number (0.54), thousand seed weight (0.48), seed width (0.49), and seed yield (0.54). EGBLUP had the best predictive ability for days to flowering (0.19), branch number (0.20), and RF had the best predictive ability for seed area (0.56), and seed length (0.56). The standard GBLUP model did not outperform for any traits, but its magnitude was very close to that of the corresponding best predictive ability (Fig. 3 and Supplementary Fig. S7). We also calculated the correlation between GEBVs estimated by 13 (SVM was not considered) different models to observe the consistency between them (Fig. 4 and Supplementary Fig. S8). In the case of days to flowering and branch number, model LASSO and EN had less correlation (r ≤ 0.85) with all other models. Other than this, for all traits, all models showed a very good correlation (r ≥ 0.92) while many pairwise correlations were close to one. These results indicate that all models excluding LASSO and EN yielded highly consistent GEBVs across traits.
Predictive ability considering various marker subsets
The predictive ability for different traits varies according to marker density. For all traits, the predictive ability yielded by marker subset (5362 markers), based on LD decay distance, was almost similar to that by the whole set of markers, with the exception of days to flowering. The predictive ability was least for 20 randomly selected markers, and it increased to a steady plateau above 3000 randomly selected markers for all traits (Table 3).
| Marker selection criteria | No. of markers | Predictive ability | |||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| DF | PH | TL | BN | BollN | TSW | SA | SW | SL | Yield | ||
| All marker | 26171 | 0.19 | 0.48 | 0.56 | 0.20 | 0.54 | 0.48 | 0.56 | 0.49 | 0.56 | 0.54 |
| All marker* | 26171 | 0.23 | 0.52 | 0.60 | 0.38 | 0.63 | 0.50 | 0.57 | 0.53 | 0.59 | 0.58 |
| LD pruning | 5362 | 0.15 | 0.48 | 0.56 | 0.21 | 0.54 | 0.48 | 0.55 | 0.49 | 0.54 | 0.54 |
| RS | 13000 | 0.19 | 0.48 | 0.56 | 0.20 | 0.54 | 0.49 | 0.55 | 0.49 | 0.55 | 0.54 |
| RS | 7000 | 0.18 | 0.48 | 0.56 | 0.20 | 0.54 | 0.48 | 0.56 | 0.48 | 0.54 | 0.54 |
| RS | 3000 | 0.17 | 0.48 | 0.55 | 0.19 | 0.54 | 0.47 | 0.55 | 0.48 | 0.55 | 0.54 |
| RS | 1000 | 0.14 | 0.49 | 0.56 | 0.20 | 0.54 | 0.48 | 0.54 | 0.48 | 0.54 | 0.54 |
| RS | 200 | 0.14 | 0.45 | 0.54 | 0.15 | 0.51 | 0.45 | 0.52 | 0.46 | 0.52 | 0.50 |
| RS | 20 | 0.04 | 0.42 | 0.42 | 0.06 | 0.38 | 0.34 | 0.42 | 0.31 | 0.27 | 0.33 |
We selected marker subsets based on GWAS and one-way ANOVA results using the whole genotype set as well as the training set only (Table 4). The marker subsets selected, based on scenario-I, yielded the lowest (0.28) predictive ability for days to flowering, branch number, and best predictive ability (0.72) for thousand seed weight. Scenario-II showed a predictive ability ranging from 0.45 (branch number) to 0.67 (seed length). In the case of scenario-III, the predictive ability for all traits significantly dropped compared to that of scenario-I & II. The predictive ability computed using all markers was better than that of scenario-III but was lower than that of scenario-I & II for all traits (Table 4). At all p-levels, the predictive ability for scenario-IV was better than that of scenario-V for days to flowering, technical length, branch number, and thousand seed weight, but opposite results were found for plant height, boll number, seed area, seed width, seed length, and seed yield. The difference between predictive ability based on scenario-IV & V was better for days to flowering, branch number, seed area, seed length and seed yield. For branch number and days to flowering, predictive ability based on scenario-IV was better than that by the whole set of markers, but for scenario-V, it was smaller than that by the whole set of markers. Except for days to flowering and branch number, for all other traits, the predictive ability at the corresponding p-level for scenario-IV & V was almost similar to that of the whole marker set (Table 4). The marker number used for predictive ability computation varied according to traits and selection scenarios (Supplementary Table S3).
| Marker selection criteria | P-value | Predictive ability | |||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| DF | PH | TL | BN | BollN | TSW | SA | SW | SL | Yield | ||
| All marker | - | 0.19 | 0.48 | 0.56 | 0.20 | 0.54 | 0.48 | 0.56 | 0.49 | 0.56 | 0.54 |
| GWAS_WSa | 0.000103 | 0.28 | 0.58 | 0.63 | 0.28 | 0.58 | 0.72 | 0.61 | 0.52 | 0.67 | 0.61 |
| GWAS_WSb | - | 0.47 | 0.60 | 0.61 | 0.45 | 0.61 | 0.64 | 0.51 | 0.60 | 0.67 | 0.64 |
| GWAS_TSc | - | 0.02 | 0.27 | 0.17 | 0.14 | 0.39 | 0.33 | 0.41 | 0.33 | 0.40 | 0.27 |
| ANOVA_WSd | 0.001 | 0.45 | 0.49 | 0.55 | 0.47 | 0.53 | 0.51 | 0.46 | 0.52 | 0.47 | 0.53 |
| ANOVA_TSe | 0.001 | 0.05 | 0.51 | 0.52 | 0.15 | 0.56 | 0.46 | 0.53 | 0.51 | 0.51 | 0.58 |
| ANOVA_WSd | 0.01 | 0.49 | 0.48 | 0.54 | 0.45 | 0.52 | 0.49 | 0.44 | 0.49 | 0.45 | 0.52 |
| ANOVA_TSe | 0.01 | 0.11 | 0.51 | 0.53 | 0.11 | 0.56 | 0.46 | 0.53 | 0.51 | 0.51 | 0.58 |
| ANOVA_WSd | 0.05 | 0.48 | 0.47 | 0.53 | 0.44 | 0.51 | 0.48 | 0.43 | 0.48 | 0.43 | 0.51 |
| ANOVA_TSe | 0.05 | 0.14 | 0.50 | 0.53 | 0.12 | 0.55 | 0.46 | 0.53 | 0.51 | 0.51 | 0.58 |
Predictive ability considering population structure
We partitioned the whole collection into five (P1 to P5) (Fig. 1A) and four (C1 to C4) (Fig. 1B) clusters based on marker genotype and phenotypic data, respectively. To increase the level of relatedness among genotypes, genotypes belonging to cluster P3 and (P3 + P4) were discarded from the whole collection as these two clusters showed the greatest divergence from other clusters and incorporated Q-matrix in the model. In all cases, predictive ability was lower than that resulting from using the whole collection for all traits (Table 5). Predictive ability within each cluster was always lower than that of using the whole collection for all traits. The magnitude of genotype number within each cluster was not reflected by the magnitude of predictive ability i.e., clusters having more genotypes exhibited both high and low predictive ability for different traits and vice-versa (Table 5).
| Clusters | Genotype number | Predictive ability | |||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| DF | PH | TL | BN | BollN | TSW | SA | SW | SL | Yield | ||
| Whole set (WS) | 337 | 0.19 | 0.48 | 0.56 | 0.20 | 0.54 | 0.48 | 0.56 | 0.49 | 0.56 | 0.54 |
| WS-P3a | 312 | 0.19 | 0.42 | 0.52 | 0.17 | 0.52 | 0.49 | 0.55 | 0.46 | 0.55 | 0.52 |
| WS-(P3 + P4)b | 280 | 0.18 | 0.33 | 0.40 | 0.18 | 0.43 | 0.49 | 0.52 | 0.45 | 0.53 | 0.42 |
| SP+ | 337 | 0.17 | 0.36 | 0.34 | 0.07 | -0.10 | 0.40 | 0.41 | 0.36 | 0.45 | -0.03 |
| P1 | 41 | 0.01 | -0.32 | -0.37 | 0.33 | -0.16 | 0.48 | 0.42 | 0.35 | 0.38 | -0.15 |
| P2 | 83 | 0.07 | -0.01 | -0.01 | 0.09 | -0.04 | 0.04 | 0.08 | 0.06 | 0.09 | 0.17 |
| P3 | 25 | -0.27 | -0.44 | -0.47 | -0.25 | -0.44 | 0.34 | 0.49 | 0.46 | 0.30 | -0.43 |
| P4 | 32 | 0.07 | -0.28 | -0.36 | -0.19 | 0.22 | 0.12 | 0.11 | 0.18 | -0.02 | 0.0002 |
| P5 | 156 | 0.22 | 0.11 | 0.14 | 0.11 | 0.32 | 0.22 | 0.13 | 0.09 | 0.18 | 0.19 |
| C1 | 43 | -0.29 | -0.25 | -0.09 | 0.10 | -0.28 | -0.23 | 0.10 | -0.12 | 0.06 | 0.04 |
| C2 | 156 | 0.18 | -0.001 | 0.20 | 0.14 | 0.25 | 0.08 | 0.15 | 0.09 | 0.15 | 0.14 |
| C3 | 105 | 0.20 | 0.22 | 0.28 | 0.13 | 0.25 | 0.31 | 0.19 | 0.06 | 0.07 | 0.21 |
| C4 | 33 | 0.09 | 0.27 | 0.20 | -0.29 | -0.26 | 0.25 | 0.05 | 0.03 | 0.14 | -0.25 |
Stratified sampling and calculated predictive ability were done for all traits according to M-I and M-II. In the case of M-I, predictive ability based on different sample sizes was better than that by the whole set (Fig. 5A). Sample size 125 yielded the best predictive ability for branch number, boll number, thousand seed weight, seed width, and seed length. For days to flowering, plant height, and technical length, we found the best predictive ability at sample size 50, whereas it was better for seed length and seed yield at sample size 175 and 75, respectively. Likewise, M-I, as in the case of M-II, predictive ability was better based on different sample sizes than that for the whole set across traits (Fig. 5B). For days to flowering, plant height, technical length, branch number, and boll number the predictive ability reached a plateau at sample size 175, whereas, for thousand seed weight, seed area, seed width, seed length, and seed yield the same happened at sample size 200. In most cases, predictive ability based on M-I was better than that based on M-II.
Indirect predictive ability
Predictive ability of any target trait was calculated using the GEBVs of another trait (Table 6). Trait combinations with the best positive correlation resulted in better predictive ability, except for the combinations of days to flowering and plant height, days to flowering and branch number, and plant height and branch number. For instance, the predictive ability for seed yield based on GEBVs of plant height, technical length, and boll number was 0.46, 0.50, and 0.55, respectively, as seed yield showed a better correlation to plant height (0.25), technical length (0.38), and boll number (0.73). We found similar results for other traits also.
| Traits | DF | PH | TL | BN | BollN | TSW | SA | SW | SL | Yield |
|---|---|---|---|---|---|---|---|---|---|---|
| DF | 0.19 | -0.16 | -0.16 | -0.09 | -0.01 | -0.05 | -0.07 | -0.08 | -0.02 | 0.00 |
| PH | 0.36 | 0.48 | 0.55 | 0.11 | 0.43 | -0.20 | -0.27 | -0.15 | -0.30 | 0.46 |
| TL | 0.07 | 0.78 | 0.56 | 0.10 | 0.48 | -0.21 | -0.30 | -0.17 | -0.33 | 0.50 |
| BN | 0.46 | 0.69 | 0.37 | 0.20 | 0.15 | 0.21 | 0.21 | 0.26 | 0.17 | 0.15 |
| BollN | -0.01 | 0.27 | 0.31 | 0.14 | 0.54 | -0.19 | -0.31 | -0.18 | -0.32 | 0.55 |
| TSW | -0.04 | -0.21 | -0.21 | -0.12 | -0.04 | 0.48 | 0.55 | 0.48 | 0.55 | -0.22 |
| SA | -0.01 | -0.18 | -0.21 | -0.05 | -0.17 | 0.3 | 0.56 | 0.49 | 0.57 | -0.27 |
| SW | -0.05 | -0.16 | -0.10 | -0.10 | -0.06 | 0.70 | 0.84 | 0.49 | 0.53 | -0.20 |
| SL | 0.02 | -0.24 | -0.22 | -0.12 | -0.14 | 0.74 | 0.91 | 0.82 | 0.56 | -0.29 |
| Yield | -0.01 | 0.25 | 0.38 | 0.06 | 0.73 | -0.07 | -0.25 | -0.12 | -0.20 | 0.54 |
Phenotypic variability
Here, we investigated seed yield and nine other agronomic traits such as days to flowering, plant height, technical length, branch number, boll number, thousand seed weight, seed area, seed width, and seed length. These traits play an important role in flax development, adaptation, domestication, and improvement 113. All these traits had continuous variation in all environments suggesting polygenic inheritance. Correlations among polygenic traits occur due to linkage and/or pleiotropic effect 114. The better the correlations among traits the more likely it is that breeders can indirectly select one trait based on other traits with high heritability. In the current study, we found a very good positive correlation between plant height and technical length, and among seed-related traits. This finding will allow breeders to phenotype only plant height or technical length and any seed-related traits for further research among this germplasm collection, which will greatly reduce the phenotyping and analysis workload. Among the traits, boll number had the best correlation with seed yield. Previous studies 115,116 also confirmed the most direct contribution of boll or fruit number to flax seed yield. Days to flowering always had a negative correlation to seed yield and seed-related traits, which directs to assume that there is possibility of exhausting more photosynthetic carbohydrates by late flowering flax genotypes for vegetative growth rather than seed formation and development. Seed traits such as thousand seed weight, seed area, seed width, and seed length were also negatively correlated to seed yield, which is consistent with previous studies 117–119.
Population structure
Structural variation and phenotypic diversity in a collection is inevitable when a breeder deals with a germplasm collection of different origins, types, and sources. Structural variation occurs due to allelic diversity present in the collection, whereas phenotypic variation is linked to this allelic diversity as well as environmental variations. Structure presence in a population influences its conservation and utilization and affects the output of genome-wide association analysis and genomic prediction. To choose divergent genotypes as parents and to incorporate the effect of structure in subsequent analyses, we analyzed the SNP-based diversity of the present germplasm collection. We used both the Delta K approach 76 and four alternative statistics 77 to identify the optimum number of clusters. Though the Delta K approach is widely used, it gives a variable number of clusters at different combinations of burn-in lengths and MCMC lengths. We found this limitation of the Delta K approach in this study as well as in previous studies 120,121. To overcome this limitation, four alternative statistics methods 77 were used, which revealed a consistent number (4 to 5) of clusters at different combinations of burn-in lengths and MCMC lengths. Based on this output, the germplasm collection was partitioned into five clusters. Population structure is influenced by mating strategy, mutation, selection, and gene flow 122. The clear-cut separation of Asian (P3) and Turkish (P4) genotypes from others indicates that the geographic distance accelerates genetic differentiation by hindering gene flow. The presence of variable types of genotypes having mixed origin in sub-population P2 supports the hypothesis of active exchange of germplasm among European countries 123 as well as among other countries. The grouping of all NDSU-released varieties and advanced breeding lines and Canadian genotypes under the same sub-population P5 was due to shared ancestors and exchange of germplasm between the USA and Canada 124. The fiber-type genotypes were in this sub-population as they were part of the parental set used for developing advanced breeding lines. We also partitioned the germplasm collection into four groups using ten quantitative agronomic traits. However, no clear-cut phenotypic clustering pattern according to types and origins was observed; the spring-type NDSU advanced breeding lines dominate one group, whereas spring-types of other origins and fiber type cluster together in another group. The pattern of the genotypic clusters was not reflected by phenotypic grouping and vice versa (Fig. 1). The mismatch between genotypic and phenotypic clustering output was also reported in flax 125, winged yam 126, and durum wheat 127. This mismatch could be improved by incorporating more plant features i.e., traits (qualitative and quantitative) and diverse environments in further studies.
Genomic prediction
The studied germplasm collection showing considerable genetic diversity will resist genetic erosion, boost genetic gain, and serve as a source of valuable genes for further improvement. As the studied traits had continuous variations, relying on phenotypic evaluation alone can be expensive, laborious and time intensive. However, the low genotyping cost, relatively accurate genotyping and efficient computational algorithms increase the opportunities to evaluate and utilize this collection for genomic prediction, which will reduce the time and cost associated with traits evaluation 128–132.
Comparing genomic prediction models
Various genomic prediction models are available, which can capture both linear (additive) and non-linear (epistasis and dominance) effects. In this study, we used prediction ability and time requirement to compare 14 different models. Breeders can also use prediction accuracy to compare the models proposed by Ould Estaghvirou et al. (2013) 112. Although prediction accuracy across all traits was better than corresponding prediction ability values, we chose prediction ability as criteria to compare models, since there are possibilities of estimating and interpreting heritability poorly 133,134.
Among all models, SVM yielded the lowest predictive ability for all traits. A similar performance of SVM was found in wheat 34,135, but the opposite scenario in maize 136. Apart from this, we found close predictive ability values by different models across traits, similar to that reported previously 35,36,42. For example, Bari et al. (2021) 36 determined the predictive ability for six traits in a diverse pea germplasm collection using five different models. They observed almost the same predictive ability values across traits by different models. We also found no single model worked best for all traits and this output is in line with the research findings of other authors. For instance, Azodi et al. (2019) 35, investigated 12 linear and non-linear models for different traits in six species and concluded that no single model worked best for all traits across species. Yu et al. (2022) 51, Phumichai et al. (2022) 137, and Roorkiwal et al. (2016) 138 revealed the same phenomenon in rice, cassava, and chickpea, respectively.
Along with prediction ability, breeders can compare different models’ efficiency based on the rankings of GEBVs of target genotypes and run-time requirements by those models. In our study, considering the very close correlation among GEBVs and high variations in time requirement by different models across traits, we can suggest historic GBLUP as a standard method for its simplicity and speed. Models’ efficiency is also affected by the complex and unique underlying biology behind the traits 129, 139. Poor predictive ability was observed for days to flowering and branch number compared to other traits by all models. Compared to our findings, Lan et al. (2020) 62 found better predictive ability for days to maturity, but lower predictive ability for seed yield in a bi-parental flax population of 260 lines. The lower predictive ability for seed yield was also supported by You et al. (2016) 63 in three different bi-parental flax populations. In our study, the predictive ability values by different models across traits were proportionate to traits’ heritability. It indicates heritability as a reason behind this difference in predictive ability across traits by different models. Previous reports also confirmed that heritability affects genomic prediction in many crops 42,43,140. Gene action (additive, dominance, and epistasis) also affects prediction accuracy for traits 141,142. Empirical and theoretical evidence indicates that the lion share of genetic variance is additive, though gene action is not 143,144. Momen et al. (2018) 37 reported that linear or parametric, and non-linear or non-parametric models outperform for traits under additive and non-additive gene action, respectively. In our study, outperforming the linear RR model for most of the traits conferred that these traits were under additive gene action. Apart from this, linear model EGBLUP and non-linear model RF yielded the best prediction ability for days to flowering, branch number, seed area, and seed length. Epistatic gene action may shape these traits as outperforming models were supposed to capture epistatic genetic variance 37,91.
Marker density effect on predictive ability
Cost-effective next-generation sequencing techniques and availability of high-quality reference genome have enabled breeders to extract informative genetic markers in prolific numbers. Utilizing these resources and high-performance computing facilities, breeders can feed the models with a huge number of markers. Although many models have been developed to handle the over-parameterization problem (marker number > > observation number) in genomic selection, previous reports confirmed that adding more markers after a certain number did not improve prediction ability. For example, in a wheat panel of 760 lines, prediction ability reached a plateau above around 5000 randomly selected markers 34, whereas in maize natural and bi-parental populations it requires 7000 and 2000 randomly selected markers, respectively 141. We found the same trend of prediction ability in this study where a plateau was obtained around 1000 to 3000 randomly selected markers. In a recent simulation study, Chang et al. (2018) 44 achieved the same prediction accuracy by using 0.5 to 1% of all markers compared to that of the whole marker set (200,000). In our case, it happened at around 26% (7000) of total markers. This finding indicates that the prediction ability by marker subset capturing all QTL information and by the whole marker set will be the same. Our finding confirmed this, where prediction ability by a marker subset (5362) based on LD-decay distance was similar to that by the whole marker set. In this study, more markers were required to obtain maximum prediction ability, though it required only 256 markers in a wheat bi-parental population 145, and 1000–1200 markers in a soybean varietal collection of 235 individuals 146. This discrepancy in marker numbers among various research was due to the nature of the studied population and LD-decay pattern. Although fewer markers can capture all QTL information in a bi-parental and varietal collection due to slow LD-decay, more markers are required in a germplasm collection due to rapid LD-decay.
Marker-trait association effect on predictive ability
In this study, the GWAS-derived significant SNP subset yielded better predictive ability values than that by the randomly selected marker subset and, surprisingly, even better predictive ability than that by the whole marker set across all traits. This finding was in line with results obtained by other authors in flax 62, wheat 34,45, maize 46,147 and spinach 148. This overestimation only prevailed when marker selection was made on the whole population (training set + validation set), instead of the training set only, indicating a clear overfitting. That is why in practice, using the output by whole population-based significant markers is not recommended. This overfitting did not occur when marker selection was done based on a one-way ANOVA approach either using the whole collection or training set across all traits, except for days to flowering and branch number. This was due to capturing more markers (≥ 2000) by the one-way ANOVA approach, as more markers dissolve overfitting by acting like random selection. For future studies, breeders can use marker subsets based on LD-decay distance rather than random selection and GWAS-based selection, which will ensure maximum prediction ability and will reduce run time.
Population structural effect on predictive ability
Because availability of diverse genetic materials in a breeding program ensures its sustainability, breeders strive to improve predictive ability while maintaining genetic diversity. The studied genotypic collection has substantial genotypic and phenotypic diversity (Fig. 1). We corrected structural variation by discarding most divergent clusters and incorporating the population Q-matrix during analysis. In both cases, predictive ability did not improve across traits. Similar phenomena were found in wheat 149,150, pea 36,151, barley 152 and maize 47,153. This finding confirms that for quantitative traits, reducing genetic structural variation increases genetic homogeneity in the collection, but not phenotypic homogeneity. We also observed very low to negative predictive ability across traits within each genotypic and phenotypic cluster, though Haile et al. (2021) 154 found moderate to high prediction accuracies within wheat subpopulations. The reason behind this was the inconsistency between genotypic and phenotypic variance i.e., genotypes of specific genetic clusters were grouped under different phenotypic clusters and the same happened to the genotypes of specific phenotypic clusters (Fig. 1C-D). Smaller population sizes in some clusters may also contribute to this.
In this study, the reduction of predictive ability due to the correction of population structure has driven us to follow a stratified sampling approach, as previous research indicates that stratified sampling 48 or composite sampling 155 increases predictive ability. We made a stratified sample of different sizes (25 to 300) following methods M-I and M-II as described in the materials and methods section. In agreement with previous studies, stratified sampling yielded better predictive ability compared to the whole collection. For all traits, in the case of M-I, a small sample size (50 to 125) yielded better predictive ability compared to the whole collection. In the case of M-II, predictive ability increased with a larger training set size, indicating M-II behaves like random sampling. In terms of predictive ability, M-I was more productive than M-II because M-I always captures fewer genotypes than M-II. These results confirmed that small sized stratified sample is strong enough to capture diversity as well as maintain good predictive ability in a diverse germplasm collection. In addition, breeders can use genotypes as parents from a particular sample yielding the increased predictive ability to make Multi-parent Advanced Generation Intercrosses (MAGIC) populations for future studies because the relatedness among individuals of the training set and target set accelerates predictive ability 38,39,156.
Indirect genomic prediction
In this study, traits had both positive and negative correlations with each other. Breeders can utilize information on correlated traits to predict the target trait using a multi-trait genomic prediction approach. Many previous studies exhibited benefit 40,41,157,158 and no benefit 159–161 of multi-trait over single-trait prediction approach in different crops. However, one of the major limitations of the multi-trait approach is that breeders need to phenotype multiple correlated traits, which is expensive and laborious. To overcome this limitation, breeders can use an indirect genomic prediction approach i.e., prediction of genotype for target trait by using a correlated single trait. The benefit of the indirect approach accelerates if the correlated trait possesses high heritability and is easy to phenotype at crop early stages of development and vice-versa for the focal trait. We found indirect prediction abilities were generally proportional to the correlation coefficient of traits with target traits. For example, the indirect prediction abilities for seed yield by plant height, technical length and boll number increased gradually with the increment of correlation coefficient values. For all traits, indirect predictive ability based on highly correlated traits was very close to its single-trait predictive ability. A similar result was found by Fernandes et al. (2018) 161 in sorghum where they reported the indirect prediction accuracy for biomass yield by plant height was similar to its single-trait and multi-trait prediction accuracies. Our findings will help breeders to reduce workload by performing indirect selection for expensive or labor-intensive focal traits by phenotyping early expressed correlated simple-to-measure traits at an early stage of the breeding pipeline.
In this study, a rigorous investigation of various key factors affecting genomic prediction ability was conducted by comparing fourteen different models. The results indicated almost similar performance by models in terms of predictive ability, and no single model worked best across all traits. Predictive ability reaches a plateau around a certain marker density and shows similarity with the whole marker set when choosing markers covering all QTL information. In the diverse flax collection used for this study, the small sample size representing population structure was strong enough to boost predictive ability compared to the whole collection across all traits. Along with this, indirect selection for seed yield considering correlated traits also holds potential for applied breeding efforts. This research presented herein will equip the plant breeders to efficiently design various aspects of genomic prediction to gain increased selection accuracy, which will subsequently accelerate the program’s genetic gain.
Acknowledgments: We thank Mr. Kreg Kercher (Department of Plant Sciences, NDSU, Fargo, ND) for helping in field planting and data collection. We also thank Justin D Faris (Research Geneticist, USDA-ARS, Fargo, ND) for providing MARVIN seed analyzer. Technical assistance during data analyses from Md. Abdullah Al Bari (Department of Plant Sciences, NDSU, Fargo, ND), Jason D. Fiedler (Research Plant Molecular Geneticist, USDA-ARS, Fargo, ND) and Brant Bigger (Plant Physiologist, USDA-ARS, Fargo, ND) is gratefully acknowledged. This work used resources of the Center for Computationally Assisted Science and Technology (CCAST) at North Dakota State University, which were made possible in part by NSF MRI Award No. 2019077.
Author contributions: AH and MR conceived and designed the study. AH conducted the data collection, curation, and analyses. The interpretation of results was done by AH with the help of MR and JVA. AH wrote the manuscript. MR is the principal investigator of the project. All authors participated in revising and editing the manuscript and approved the final version of the manuscript.
Data Availability: All raw sequence data and variant data are available in NCBI and EVA repository. The accession id for them is PRJNA979944 (https://www.ncbi.nlm.nih.gov/sra/PRJNA979944) and PRJEB62432 (https://www.ebi.ac.uk/eva/?eva-study=PRJEB62432), respectively. The Phenotypic datasets and R scripts used and/or analyzed during the current study are available from the corresponding author upon reasonable request.
Funding: The study was funded by the U.S. Department of Agriculture - National Institute of Food and Agriculture (Hatch Project No. ND01581). The funders had no role in study design, data collection and analysis, the decision to publish, or preparation of the manuscript.
Competing Interests: On behalf of all authors, the corresponding author states that there is no conflict of interest.
Ethics approval: Not applicable.
Consent to participate: Not applicable.
Consent to publish: Not applicable.
- Muir, A. D. & Westcott, N. D. Flax: the genus Linum. (CRC press, 2003).
- Touré, A. & Xueming, X. Flaxseed Lignans: Source, Biosynthesis, Metabolism, Antioxidant Activity, Bio-Active Components, and Health Benefits. Compr Rev Food Sci Food Saf 9, 261–269 (2010). https://doi.org/10.1111/J.1541-4337.2009.00105.X
- Westcott, N. D. & Muir, A. D. Flax seed lignan in disease prevention and health promotion. Phytochemistry Reviews 2, 401–417 (2003).
- Przybylski, R. Flax oil and high linolenic oils. Bailey’s industrial oil and fat products. 2:281–301 (2005).
- Morris, D. H. The novel egg: Opportunities for flax in omega-3 egg production. Flax Council of Canada (2003). http://www.flaxcouncil.ca/english/pdf/novelegg.pdf
- Ndou, S. P., Kiarie, E., Walsh, M. C. & Nyachoti, C. M. Nutritive value of flaxseed meal fed to growing pigs. Anim Feed Sci Technol 238, 123–129 (2018).
- Cullis, C. Linum. in Wild crop relatives: Genomic and breeding resources (ed. Kole, C.) 177–189 (Springer, 2011).
- FAOSTAT. Food and agriculture organization of the United Nations. Statistical database (2022).
- USDA-NASS. United States Census of Agriculture (2022).
- Berglund, D. R. & Zollinger, R. K. Flax production in North Dakota. North Dakota Agricultural Experimental Station, Extension Service North Dakota, Fargo, North Dakota, USA. Bull. A-1038. 12 p. (2002).
- Duguid, S., Lafond, G., McAndrew, D. W., Rashid, K. Y. & Ulrich, A. Growing Flax: Production, Management & Diagnostic Guide. Winnipeg, MB: Flax Council of Canada (2007).
- Xu, Y. & Crouch, J. H. Marker-Assisted Selection in Plant Breeding: From Publications to Practice. Crop Sci 48, 391–407 (2008).
- Nihad, S. A. I. et al. Linkage of SSR markers with rice blast resistance and development of partial resistant advanced lines of rice (Oryza sativa) through marker-assisted selection. Physiology and Molecular Biology of Plants 1–17 (2022) doi:10.1007/S12298-022-01141-3.
- Sun, L. et al. Robust identification of low-Cd rice varieties by boosting the genotypic effect of grain Cd accumulation in combination with marker-assisted selection. J Hazard Mater 424, 127703 (2022).
- Alsaleh, A. et al. Marker-assisted selection and validation of DNA markers associated with cadmium content in durum wheat germplasm. Crop Pasture Sci (2022) doi:10.1071/CP21484.
- Soriano, M. et al. Identification and Characterisation of Stripe Rust Resistance Genes Yr66 and Yr67 in Wheat Cultivar VL Gehun 892. Agronomy 12, 318 (2022).
- Yadav, P. S. et al. Enhanced resistance in wheat against stem rust achieved by marker assisted backcrossing involving three independent Sr genes. Curr Plant Biol 2, 25–33 (2015).
- Yang, R., Yan, Z., Wang, Q., Li, X. & Feng, F. Marker-assisted backcrossing of lcyE for enhancement of proA in sweet corn. Euphytica 214, 1–12 (2018).
- Hao, X., Li, X., Yang, X. & Li, J. Transferring a major QTL for oil content using marker-assisted backcrossing into an elite hybrid to increase the oil content in maize. Molecular Breeding 34, 739–748 (2014).
- Yathish, K. R. et al. Introgression of the Low Phytic Acid Locus (lpa2) Into Elite Maize (Zea Mays L.) Inbreds Through Marker-Assisted Backcross Breeding (MABB). Euphytica 218, 127 (2022). doi:10.21203/rs.3.rs-1293507/v1.
- Becker, H. C. & Bernardo, R. A model for marker-assisted selection among single crosses with multiple genetic markers. Theoretical and Applied Genetics 97, 473–478 (1998).
- Bernardo, R. BREEDING for QUANTITATIVE TRAITS in PLANTS Third Edition. Stemma Press, Woodbury, MN (2020).
- Meuwissen, T. H. E., Hayes, B. J. & Goddard, M. E. Prediction of Total Genetic Value Using Genome-Wide Dense Marker Maps. Genetics 157, 1819–1829 (2001).
- Lorenz, A. J. et al. Genomic Selection in Plant Breeding: Knowledge and Prospects. Advances in Agronomy 110, 77–123 (2011).
- Schaeffer, L. R. Strategy for applying genome-wide selection in dairy cattle. Journal of Animal Breeding and Genetics 123, 218–223 (2006).
- Endelman, J. B. Ridge Regression and Other Kernels for Genomic Selection with R Package rrBLUP. Plant Genome 4, 250–255 (2011).
- Long, N., Gianola, D., Rosa, G. J. M. & Weigel, K. A. Application of support vector regression to genome-assisted prediction of quantitative traits. Theoretical and Applied Genetics 123, 1065–1074 (2011).
- de Los Campos, G. et al. Predicting Quantitative Traits with Regression Models for Dense Molecular Markers and Pedigree. Genetics 182, 375–385 (2009).
- Park, T. & Casella, G. The Bayesian Lasso. Journal of the American Statistical Association 103, 681–686 (2012). https://doi.org/10.1198/016214508000000337
- Crossa, J. et al. Prediction of genetic values of quantitative traits in plant breeding using pedigree and molecular markers. Genetics 186, 713–724 (2010).
- Ober, U. et al. Predicting Genetic Values: A Kernel-Based Best Linear Unbiased Prediction With Genomic Data. Genetics 188, 695–708 (2011).
- Long, N., Gianola, D., Rosa, G. J. M. & Weigel, K. A. Marker-assisted prediction of non-additive genetic values. Genetica 139, 843–854 (2011).
- de Los Campos, G., Gianola, D., Rosa, G. J. M., Weigel, K. A. & Crossa, J. Semi-parametric genomic-enabled prediction of genetic values using reproducing kernel Hilbert spaces methods. Genet Res (Camb) 92, 295–308 (2010).
- Charmet, G., Tran, L. G., Auzanneau, J., Rincent, R. & Bouchet, S. BWGS: A R package for genomic selection and its application to a wheat breeding programme. PLoS One 15, e0222733 (2020).
- Azodi, C. B. et al. Benchmarking Parametric and Machine Learning Models for Genomic Prediction of Complex Traits. G3 Genes|Genomes|Genetics 9, 3691–3702 (2019).
- Bari, M. A. al et al. Harnessing Genetic Diversity in the USDA Pea Germplasm Collection Through Genomic Prediction. Front Genet 12, 2273 (2021).
- Momen, M. et al. Predictive ability of genome-assisted statistical models under various forms of gene action. Scientific Reports 8:1 8, 1–11 (2018).
- Riedelsheimer, C., Brotman, Y., Méret, M., Melchinger, A. E. & Willmitzer, L. The maize leaf lipidome shows multilevel genetic control and high predictive value for agronomic traits. Scientific Reports 3, 1–7 (2013).
- Rutkoski, J. et al. Efficient Use of Historical Data for Genomic Selection: A Case Study of Stem Rust Resistance in Wheat. The Plant Genome 8, 1 (2015). 10.3835/plantgenome2014.09.0046
- Bhatta, M. et al. Multi-trait Genomic Prediction Model Increased the Predictive Ability for Agronomic and Malting Quality Traits in Barley (Hordeum vulgare L.). G3 Genes|Genomes|Genetics 10, 1113–1124 (2020).
- Velazco, J. G. et al. Genomic Prediction of Grain Yield and Drought-Adaptation Capacity in Sorghum Is Enhanced by Multi-Trait Analysis. Front Plant Sci 10, 997 (2019).
- Spindel, J. et al. Genomic Selection and Association Mapping in Rice (Oryza sativa): Effect of Trait Genetic Architecture, Training Population Composition, Marker Number and Statistical Model on Accuracy of Rice Genomic Selection in Elite, Tropical Rice Breeding Lines. PLoS Genet 11, e1004982 (2015).
- Zhang, A. et al. Effect of trait heritability, training population size and marker density on genomic prediction accuracy estimation in 22 bi-parental tropical maize populations. Front Plant Sci 8, 1916 (2017).
- Chang, L. Y., Toghiani, S., Ling, A., Aggrey, S. E. & Rekaya, R. High density marker panels, SNPs prioritizing and accuracy of genomic selection. BMC Genet 19, 1–10 (2018).
- Sehgal, D. et al. Incorporating Genome-Wide Association Mapping Results into Genomic Prediction Models for Grain Yield and Yield Stability in CIMMYT Spring Bread Wheat. Front Plant Sci 11, 197 (2020).
- Rice, B. & Lipka, A. E. Evaluation of RR-BLUP Genomic Selection Models that Incorporate Peak Genome-Wide Association Study Signals in Maize and Sorghum. Plant Genome 12, 180052 (2019).
- Guo, Z. et al. The impact of population structure on genomic prediction in stratified populations. Theor Appl Genet 127, 749–762 (2014).
- Isidro, J. et al. Training set optimization under population structure in genomic selection. Theor Appl Genet 128, 145–158 (2015).
- Frouin, J., Labeyrie, A., Boisnard, A., Sacchi, G. A. & Ahmadi, N. Genomic prediction offers the most effective marker assisted breeding approach for ability to prevent arsenic accumulation in rice grains. PLoS One 14, e0217516 (2019).
- Monteverde, E. et al. Integrating Molecular Markers and Environmental Covariates To Interpret Genotype by Environment Interaction in Rice (Oryza sativa L.) Grown in Subtropical Areas. G3 Genes|Genomes|Genetics 9, 1519–1531 (2019).
- Yu, P. et al. Genome-wide association study and genomic prediction for yield and grain quality traits of hybrid rice. Molecular Breeding 42, 1–12 (2022).
- Huang, M. et al. Use of genomic selection in breeding rice (Oryza sativa L.) for resistance to rice blast (Magnaporthe oryzae). Molecular Breeding 39, 1–16 (2019).
- Ben-Sadoun, S. et al. Economical optimization of a breeding scheme by selective phenotyping of the calibration set in a multi-trait context: application to bread making quality. Theoretical and Applied Genetics 133, 2197–2212 (2020).
- Merrick, L. F., Herr, A. W., Sandhu, K. S., Lozada, D. N. & Carter, A. H. Utilizing Genomic Selection for Wheat Population Development and Improvement. Agronomy 12, 522 (2022).
- Song, J. et al. Practical application of genomic selection in a doubled-haploid winter wheat breeding program. Molecular Breeding 37, 1–15 (2017).
- Hu, X. et al. Effectiveness of Genomic Selection by Response to Selection for Winter Wheat Variety Improvement. Plant Genome 12, 180090 (2019).
- Robert, P. et al. Phenomic selection in wheat breeding: identification and optimisation of factors influencing prediction accuracy and comparison to genomic selection. Theoretical and Applied Genetics 135, 895–914 (2022).
- Cerrudo, D. et al. Genomic selection outperforms marker assisted selection for grain yield and physiological traits in a maize doubled haploid population across water treatments. Front Plant Sci 9, 366 (2018).
- Zhang, X. et al. Rapid Cycling Genomic Selection in a Multiparental Tropical Maize Population. G3 Genes|Genomes|Genetics 7, 2315–2326 (2017).
- Fristche-Neto, R., Akdemir, D. & Jannink, J. L. Accuracy of genomic selection to predict maize single-crosses obtained through different mating designs. Theoretical and Applied Genetics 131, 1153–1162 (2018).
- Mastrodomenico, A. T., Bohn, M. O., Lipka, A. E. & Below, F. E. Genomic Selection Using Maize Ex-Plant Variety Protection Germplasm for the Prediction of Nitrogen-Use Traits. Crop Sci 59, 212–220 (2019).
- Lan, S. et al. Genomic Prediction Accuracy of Seven Breeding Selection Traits Improved by QTL Identification in Flax. International Journal of Molecular Sciences 21, 1577 (2020).
- You, F. M., Booker, H. M., Duguid, S. D., Jia, G. & Cloutier, S. Accuracy of genomic selection in biparental populations of flax (Linum usitatissimum L.). Crop J 4, 290–303 (2016).
- He, L. et al. Evaluation of Genomic Prediction for Pasmo Resistance in Flax. International Journal of Molecular Sciences 20, 359 (2019).
- Elshire, R. J. et al. A Robust, Simple Genotyping-by-Sequencing (GBS) Approach for High Diversity Species. PLoS One 6, e19379 (2011).
- Glaubitz, J. C. et al. TASSEL-GBS: a high capacity genotyping by sequencing analysis pipeline. PLoS One 9, e90346 (2014).
- You, F. M. et al. Chromosome-scale pseudomolecules refined by optical, physical and genetic maps in flax. The Plant Journal 95, 371–384 (2018).
- Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. Nat Methods 9, 357 (2012).
- Danecek, P. et al. The variant call format and VCFtools. Bioinformatics 27, 2156–2158 (2011).
- Federer, W. T. & Crossa, J. I.4 Screening Experimental Designs for Quantitative Trait Loci, Association Mapping, Genotype-by Environment Interaction, and Other Investigations. Front Physiol 3, (2012).
- Nôžková, J. et al. Descriptor List for Flax–Linum usitatissimum L. Nitra: SPU (2011).
- Cullis, B. R., Smith, A. B. & Coombes, N. E. On the design of early generation variety trials with correlated data. Journal of Agricultural, Biological, and Environmental Statistics 11, 381–393 (2006).
- Covarrubias-Pazaran, G. Genome-Assisted Prediction of Quantitative Traits Using the R Package sommer. PLoS One 11, e0156744 (2016).
- Taiyun Wei, M. et al. Package ‘corrplot’ Title Visualization of a Correlation Matrix. (2017).
- Pritchard, J. K., Stephens, M. & Donnelly, P. Inference of Population Structure Using Multilocus Genotype Data. Genetics 155, 945–959 (2000).
- Evanno, G., Regnaut, S. & Goudet, J. Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol Ecol 14, 2611–2620 (2005).
- Puechmaille, S. J. The program structure does not reliably recover the correct population structure when sampling is uneven: subsampling and new estimators alleviate the problem. Mol Ecol Resour 16, 608–627 (2016).
- Li, Y. L. & Liu, J. X. StructureSelector: A web-based software to select and visualize the optimal number of clusters using multiple methods. Mol Ecol Resour 18, 176–177 (2018).
- Jakobsson, M. & Rosenberg, N. A. CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics 23, 1801–1806 (2007).
- Bradbury, P. J. et al. TASSEL: software for association mapping of complex traits in diverse samples. Bioinformatics 23, 2633–2635 (2007).
- Jakobsson, M., Edge, M. D. & Rosenberg, N. A. The relationship between FST and the frequency of the most frequent allele. Genetics (2013) doi:10.1534/genetics.112.144758.
- Excoffier, L. & Lischer, H. E. L. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Resour 10, 564–567 (2010).
- Scrucca, L., Fop, M., Murphy, T. B. & Raftery, A. E. mclust 5: Clustering, Classification and Density Estimation Using Gaussian Finite Mixture Models. R J 8, 289 (2016).
- Villanueva, R. A. M. & Chen, Z. J. ggplot2: Elegant Graphics for Data Analysis (2nd ed.). https://doi.org/10.1080/15366367.2019.1565254 17, 160–167 (2019).
- Zhang, C., Dong, S.-S., Xu, J.-Y., He, W.-M. & Yang, T.-L. PopLDdecay: a fast and effective tool for linkage disequilibrium decay analysis based on variant call format files. Bioinformatics 35, 1786–1788 (2018).
- Hill, W. G. & Weir, B. S. Variances and covariances of squared linkage disequilibria in finite populations. Theor Popul Biol 33, 54–78 (1988).
- Desta, Z. A. & Ortiz, R. Genomic selection: Genome-wide prediction in plant improvement. Trends Plant Sci 19, 592–601 (2014).
- Piepho, H. P. Ridge Regression and Extensions for Genome wide Selection in Maize. Crop Sci 49, 1165–1176 (2009).
- Habier, D., Fernando, R. L. & Garrick, D. J. Genomic BLUP decoded: A look into the black box of genomic prediction. Genetics 194, 597–607 (2013).
- VanRaden, P. M. Efficient Methods to Compute Genomic Predictions. J Dairy Sci 91, 4414–4423 (2008).
- Jiang, Y. & Reif, J. C. Modeling epistasis in genomic selection. Genetics 201, 759–768 (2015).
- Friedman, J., Hastie, T. & Tibshirani, R. Regularization Paths for Generalized Linear Models via Coordinate Descent. J Stat Softw 33, 1 (2010).
- Usai, M. G., Goddard, M. E. & Hayes, B. J. LASSO with cross-validation for genomic selection. Genet Res (Camb) 91, 427–436 (2009).
- Zou, H. & Hastie, T. Regularization and variable selection via the elastic net. Journal of the Royal Statistical Society Series B: Statistical Methodology 67,2 301–320 (2005).
- de los Campos, G., Hickey, J. M., Pong-Wong, R., Daetwyler, H. D. & Calus, M. P. L. Whole-genome regression and prediction methods applied to plant and animal breeding. Genetics 193, 327–345 (2013).
- Habier, D., Fernando, R. L., Kizilkaya, K. & Garrick, D. J. Extension of the bayesian alphabet for genomic selection. BMC Bioinformatics 12, 1–12 (2011).
- Pérez, P. & de Los Campos, G. Genome-wide regression and prediction with the BGLR statistical package. Genetics 198, 483–495 (2014).
- Gianola, D. & van Kaam, J. B. C. H. M. Reproducing kernel Hilbert spaces regression methods for genomic assisted prediction of quantitative traits. Genetics 178, 2289–2303 (2008).
- Breiman, L. Random Forests. Machine Learning 45:1 45, 5–32 (2001).
- Maenhout, S., de Baets, B., Haesaert, G. & van Bockstaele, E. Support vector machine regression for the prediction of maize hybrid performance. Theoretical and Applied Genetics 115, 1003–1013 (2007).
- González-Recio, O., Rosa, G. J. M. & Gianola, D. Machine learning methods and predictive ability metrics for genome-wide prediction of complex traits. Livest Sci 166, 217–231 (2014).
- Price, A. L. et al. Principal components analysis corrects for stratification in genome-wide association studies. Nat Genet 38, 904–909 (2006).
- Yu, J. et al. A unified mixed-model method for association mapping that accounts for multiple levels of relatedness. Nat Genet 38, 203–208 (2006).
- Zhang, Z. et al. Mixed linear model approach adapted for genome-wide association studies. Nature Genetics 42, 355–360 (2010).
- Li, M. et al. Enrichment of statistical power for genome-wide association studies. BMC Biol 12, 1–10 (2014).
- Wang, Q., Tian, F., Pan, Y., Buckler, E. S. & Zhang, Z. A SUPER powerful method for genome wide association study. PLoS One 9, (2014).
- Segura, V. et al. An efficient multi-locus mixed-model approach for genome-wide association studies in structured populations. Nat Genet 44, 825–830 (2012).
- Liu, X., Huang, M., Fan, B., Buckler, E. S. & Zhang, Z. Iterative Usage of Fixed and Random Effect Models for Powerful and Efficient Genome-Wide Association Studies. PLoS Genet 12, (2016).
- Huang, M., Liu, X., Zhou, Y., Summers, R. M. & Zhang, Z. BLINK: a package for the next level of genome-wide association studies with both individuals and markers in the millions. Gigascience 8, 1–12 (2019).
- Wang, J. & Zhang, Z. GAPIT Version 3: Boosting Power and Accuracy for Genomic Association and Prediction. Genomics Proteomics Bioinformatics 19, 629–640 (2021).
- Li, J. & Ji, L. Adjusting multiple testing in multilocus analyses using the eigenvalues of a correlation matrix. Heredity 2005 95:3 95, 221–227 (2005).
- Ould Estaghvirou, S. B. et al. Evaluation of approaches for estimating the accuracy of genomic prediction in plant breeding. BMC Genomics 14, 1–21 (2013).
- Zhang, J. et al. Genomic Comparison and Population Diversity Analysis Provide Insights into the Domestication and Improvement of Flax. iScience 23, 100967 (2020).
- Saltz, J. B., Hessel, F. C. & Kelly, M. W. Trait Correlations in the Genomics Era. Trends Ecol Evol 32, 279–290 (2017).
- Çopur, O. & Demirel, U. Determination of Correlation and Path Analysis among Yield Components and Seed Yield in Oil Flax Varieties (Linum usitatissimum L.). Journal of Biological Sciences (2006) doi:10.3923/jbs.2006.738.743.
- Bibi, T., Mahmood, T., Mirza, Y. & Mahmood, T. Correlation studies of some yield related traits in linseed (Linum usitatissimum L.). J. Agric. Res 51 (2013).
- Worku, N., Heslop-Harrison, J. S. & Adugna, W. Diversity in 198 Ethiopian linseed (Linum usitatissimum) accessions based on morphological characterization and seed oil characteristics. Genet Resour Crop Evol 62, 1037–1053 (2015).
- Soto-Cerda, B. J. et al. Genomic regions underlying agronomic traits in linseed (Linum usitatissimum L.) as revealed by association mapping. J Integr Plant Biol 56, 75–87 (2014).
- Soto-Cerda, B. J. et al. Assessing the agronomic potential of linseed genotypes by multivariate analyses and association mapping of agronomic traits. Euphytica 196, 35–49 (2014).
- Yu, Z., Fredua-Agyeman, R., Hwang, S.-F. & Strelkov, S. E. Molecular genetic diversity and population structure analyses of rutabaga accessions from Nordic countries as revealed by single nucleotide polymorphism markers. BMC Genomics 22, 1–13 (2021).
- Rahman, M., Hoque, A. & Roy, J. Linkage disequilibrium and population structure in a core collection of Brassica napus (L.). PLoS ONE 17(3): e0250310. (2022) doi:10.1371/journal.pone.0250310.
- Schaal, B. A., Hayworth, D. A., Olsen, K. M., Rauscher, J. T. & Smith, W. A. Phylogeographic studies in plants: problems and prospects. Mol Ecol 7, 465–474 (1998).
- Maggioni, L. Flax Genetic Resources in Europe: Ad Hoc Meeting, 7-8 December 2001, Prague, Czech Republic. (Bioversity International, 2002).
- Fu, Y.-B., Rowland, G. G., Duguid, S. D. & Richards, K. W. RAPD analysis of 54 North American flax cultivars. Crop Sci 43, 1510–1515 (2003).
- Choudhary, S. B. et al. Genetic diversity spectrum and marker trait association for agronomic traits in global accessions of Linum usitatissimum L. Ind Crops Prod 108, 604–615 (2017).
- Agre, P. et al. Phenotypic and molecular assessment of genetic structure and diversity in a panel of winged yam (Dioscorea alata) clones and cultivars. Scientific Reports 9, 1–11 (2019).
- Royo, C. et al. Understanding the relationships between genetic and phenotypic structures of a collection of elite durum wheat accessions. Field Crops Res 119, 91–105 (2010).
- Mascher, M. et al. Genebank genomics bridges the gap between the conservation of crop diversity and plant breeding. Nature Genetics 51, 1076–1081 (2019).
- Yu, X. et al. Genomic prediction of maize microphenotypes provides insights for optimizing selection and mining diversity. Plant Biotechnol J 18, 2456–2465 (2020).
- Yu, X. et al. Genomic prediction contributing to a promising global strategy to turbocharge gene banks. Nature Plants 2, 1–7 (2016).
- Li, H., Rasheed, A., Hickey, L. T. & He, Z. Fast-Forwarding Genetic Gain. Trends Plant Sci 23, 184–186 (2018).
- Crossa, J. et al. Genomic Selection in Plant Breeding: Methods, Models, and Perspectives. Trends Plant Sci 22, 961–975 (2017).
- Piepho, H. P. & Möhring, J. Computing Heritability and Selection Response from Unbalanced Plant Breeding Trials. Genetics 177, 1881–1888 (2007).
- Dudley, J. W. & Moll, R. H. Interpretation and Use of Estimates of Heritability and Genetic Variances in Plant Breeding1. Crop Sci 9, 257–262 (1969).
- Ornella, L. et al. Genomic Prediction of Genetic Values for Resistance to Wheat Rusts. Plant Genome 5, (2012).
- Zhao, W. et al. Applications of Support Vector Machine in Genomic Prediction in Pig and Maize Populations. Front Genet 11, 1537 (2020).
- Phumichai, C. et al. Genome-wide association mapping and genomic prediction of yield-related traits and starch pasting properties in cassava. Theoretical and Applied Genetics 135, 145–171 (2022).
- Roorkiwal, M. et al. Genome-enabled prediction models for yield related traits in chickpea. Front Plant Sci 7, 1666 (2016).
- Valluru, R. et al. Deleterious Mutation Burden and Its Association with Complex Traits in Sorghum (Sorghum bicolor). Genetics 211, 1075–1087 (2019).
- Hayes, B. J., Visscher, P. M. & Goddard, M. E. Increased accuracy of artificial selection by using the realized relationship matrix. Genet Res (Camb) 91, 47–60 (2009).
- Liu, X. et al. Factors affecting genomic selection revealed by empirical evidence in maize. Crop J 6, 341–352 (2018).
- Raffo, M. A. et al. Improvement of genomic prediction in advanced wheat breeding lines by including additive-by-additive epistasis. Theoretical and Applied Genetics 135, 965–978 (2022).
- Mäki-Tanila, A. & Hill, W. G. Influence of Gene Interaction on Complex Trait Variation with Multilocus Models. Genetics 198, 355–367 (2014).
- Hill, W. G., Goddard, M. E. & Visscher, P. M. Data and Theory Point to Mainly Additive Genetic Variance for Complex Traits. PLoS Genet 4, e1000008 (2008).
- Heffner, E. L. et al. Genomic Selection Accuracy using Multifamily Prediction Models in a Wheat Breeding Program. Plant Genome 4, 65–75 (2011).
- Ma, Y. et al. Potential of marker selection to increase prediction accuracy of genomic selection in soybean (Glycine max L.). Molecular Breeding 36, 1–10 (2016).
- Bian, Y. & Holland, J. B. Enhancing genomic prediction with genome-wide association studies in multiparental maize populations. Heredity 118, 585–593 (2017).
- Shi, A. et al. Genome-wide association study and genomic prediction of white rust resistance in USDA GRIN spinach germplasm. Hortic Res 9, (2022).
- Crossa, J. et al. Genomic Prediction of Gene Bank Wheat Landraces. G3 Genes|Genomes|Genetics 6, 1819–1834 (2016).
- Norman, A., Taylor, J., Edwards, J. & Kuchel, H. Optimising Genomic Selection in Wheat: Effect of Marker Density, Population Size and Population Structure on Prediction Accuracy. G3 Genes|Genomes|Genetics 8, 2889–2899 (2018).
- Burstin, J. et al. Genetic diversity and trait genomic prediction in a pea diversity panel. BMC Genomics 16, 1–17 (2015).
- Thorwarth, P. et al. Genomic prediction ability for yield-related traits in German winter barley elite material. Theoretical and Applied Genetics 130, 1669–1683 (2017).
- Lyra, D. H. et al. Controlling population structure in the genomic prediction of tropical maize hybrids. Molecular Breeding 38, 1–17 (2018).
- Haile, T. A. et al. Genomic prediction of agronomic traits in wheat using different models and cross-validation designs. Theoretical and Applied Genetics 134, 381–398 (2021).
- He, S. et al. Genomic prediction using composite training sets is an effective method for exploiting germplasm conserved in rice gene banks. Crop J (2022) doi:10.1016/J.CJ.2021.11.011.
- Lorenz, A. & Smith, K. P. Adding Genetically Distant Individuals to Training Populations Reduces Genomic Prediction Accuracy in Barley. Crop Sci 55, 2657–2667 (2015).
- Schulthess, A. W., Zhao, Y., Longin, C. F. H. & Reif, J. C. Advantages and limitations of multiple-trait genomic prediction for Fusarium head blight severity in hybrid wheat (Triticum aestivum L.). Theoretical and Applied Genetics 131, 685–701 (2018).
- Lyra, D. H. et al. Multi-trait genomic prediction for nitrogen response indices in tropical maize hybrids. Molecular Breeding 37, 1–14 (2017).
- dos Santos, J. P. R., de Castro Vasconcellos, R. C., Pires, L. P. M., Balestre, M. & von Pinho, R. G. Inclusion of Dominance Effects in the Multivariate GBLUP Model. PLoS One 11, e0152045 (2016).
- Schulthess, A. W. et al. Multiple-trait- and selection indices-genomic predictions for grain yield and protein content in rye for feeding purposes. Theoretical and Applied Genetics 129, 273–287 (2016).
- Fernandes, S. B., Dias, K. O. G., Ferreira, D. F. & Brown, P. J. Efficiency of multi-trait, indirect, and trait-assisted genomic selection for improvement of biomass sorghum. Theoretical and Applied Genetics 131, 747–755 (2018).
No competing interests reported.
- SupplementaryFigureS1.jpg
- SupplementaryFigureS2.Phenotypicvariation.tiff
- SupplementaryFigureS3.Phenotypiccorrelation.tiff
- SupplementaryFigureS4.Correlationamongenvironments.tiff
- SupplementaryFigureS5.Clustering.png
- SupplementaryFigureS6.tiff
- SupplementaryFigureS7.BoxplotAllmarkersAlltraits.tiff
- SupplementaryFigureS8.CorrelationanalysisGEBVstitle.tiff
- SupplementaryTableS1.xlsx
- SupplementaryTableS2.xlsx
- SupplementaryTableS3.xlsx
- SupplementaryInformations.docx
