Elucidating SNP-based genetic diversity and population structure of advanced breeding lines of bread wheat (Triticum aestivum L.)

Plant genetic material

One hundred forty-one advanced breeding lines (ABLs) developed at CIMMYT, Mexico using adapted cultivars from various worldwide breeding programs were used in studying genetic diversity. These lines have been selected from a more extensive set of advanced breeding lines sent to South Asia as part of the CIMMYT wheat breeding program to develop high-yielding varieties suitable for the region. The lines were introduced in 2016 as the 4th cohort of the South Asia Bread Wheat Genomic Prediction Yield Trial (SABWGPYT). The lines were developed by making crosses with diverse and high yielding parental lines as part of the CIMMYT Global Wheat Program (GWP). A set of approx. thirty-eight thousand lines were planted in small plots in Obregon during the 2015-16 crop season. The selection was made based on yield and agronomic traits. Then lines from different genetic backgrounds were carefully selected as candidate varieties suitable for further testing in South Asia. The pedigree details data of those 141 advanced breeding lines are given in Table S1.

Genotyping-by-sequencing and SNP filtering

Genomic DNA was extracted using the modified CTAB method from the fresh leaves of wheat seedlings (Dreisigacker et al., 2016). Genotyping-by-sequencing was performed in Illumina HiSeq 2500 using protocol from Poland et al. (2012). SNP calling was performed using TASSEL v5.2.6 (Bradbury et al., 2007) using the TASSEL-GBSv2 pipeline and aligned to the reference Chinese Spring Wheat Assembly (RefSeq v1.0). Beagle v4.1 with default settings was used to impute missing data. After filter criteria quality control (sample call rate >0.8, MAF ≥ 0.05, SNP call rate >0.7), 14,563 polymorphic SNPs and 141 genotypes were selected for further analysis.

Genetic analysis of SNPs

The genotypic summary of the 14,563 SNPs was obtained using the “geno summary” function of TASSEL v5.2.6 (Bradbury et al., 2007). Chromosome-wise genomic SNP distribution, minor allele frequency (MAF), observed heterozygosity, and the polymorphism information content (PIC) were performed using GBS-based SNP markers. The PIC value of each biallelic SNP marker was computed using the following method derived from (Botstein et al., 1980): ${\displaystyle PIC=1-\left(MA{F}^2+{\left(1-MAF\right)}^2\right)-\left(2MA{F}^2{\left(1-MAF\right)}^2\right)}$

Population structure analysis

The population structure of the advanced bread wheat breeding lines was inferred using the Bayesian method implemented in STRUCTURE 2.3.4 (Pritchard, Stephens & Donnelly, 2000). Ten individualistic evaluations of every k were used. The STRUCTURE was placed on 100,000 burn-in, subsequently 100,000 Markov chain Monte Carlo (MCMC) replications. The best k for the present population was determined using the evanno algorithm (Evanno, Regnaut & Goudet, 2005) implemented in Pophelper v2.3.1 (Francis, 2017). Optimal K/sub-population was identified by the delta k-value peaks across various k values. Principal component analysis (PCA) based on covariates was performed of the SNP data in Tassel 5.2.6 (Bradbury et al., 2007) and plotted using Plotly v4.9.3 (Sievert, 2020). The PCAs were used to construct dendrograms using r package ape v5.4-1 (Paradis & Schliep, 2019) using the complete linkage method for hierarchical clustering.

Genetic diversity analysis and analysis of molecular variance (AMOVA)

The selected K value from structure output was used to subdivide the advanced wheat lines into sub-populations and was utilized for AMOVA. AMOVA and genetic diversity indices were performed for individual sub-population using GeneAlEx v6.41 (Peakall & Smouse, 2006). The percentage of molecular variance among and within subgroups was calculated from AMOVA. The haploid number of migrants (Nm) was calculated using the among-population variance (Va) and within-population variance (Vw) using the formula ${\displaystyle N_m=\left[\left\{1/\left(\frac{V_a}{V_a+V_w}\right)\right\}-1\right]∕2}$

The calculation of Shannon’s Information Index (I), effective alleles (Ne), different alleles (Na), number of loci with private alleles, unbiased diversity (uh), and the haploid genetic diversity index (h) were also performed.

Na was calculated by direct count of alleles across subpopulations per loci and averaged by the arithmetic mean across loci per sub-population. Ne was calculated using expected heterozygosity by formula ${\displaystyle Ne=\frac{1}{1-He}\text{where}He=1-\sum p^2}$

Here p is the frequency of the allele. I was calculated (per locus and averaged across the number of loci) using the formula. ${\displaystyle I=-\sum p_{i}ln{p}_i}$

Where ln is the natural logarithm of p_i, i.e., frequency of ith allele, private alleles are the alleles unique to the subpopulation. The haploid genetic diversity index (h) provides the probability that the two individuals would be different and was calculated using the formula. ${\displaystyle h=1-\sum \underset{i}{\overset{2}{p}}}$

The unbiased diversity (uh) was calculated using the allele frequency and sample size (n) with a formula. ${\displaystyle uh=\frac{n}{n-1}\left(1-\sum \underset{i}{\overset{2}{p}}\right)}$

SNPs distribution on the wheat genome

A total of 14,563 bi-allelic SNPs were equitably distributed across the three genomes. The SNP HapMap and SNP positions datasets correspond to the reference Chinese Spring Wheat Assembly (RefSeq v1.0) used for present analysis have been provided through the link (https://doi.org/10.6084/m9.figshare.14273339.v2.). The Chinese Spring Wheat Assembly (RefSeq v1.0) reference genome assembly used in the present study is a widely used reference genome in hexaploid wheat. (https://urgi.versailles.inra.fr/download/iwgsc/IWGSC_RefSeq_Assemblies/v1.0/ and https://urgi.versailles.inra.fr/download/iwgsc/IWGSC_RefSeq_Annotations/v1.0/). The B-genome had a maximum of 7378 SNPs (50.66%), subsequently the A-genome with 5921 SNPs (40.66%) and the D-genome with a minimum number of SNPs, i.e., 1264 (8.68%) (Table 1, Fig. 1A). With 14,563 markers across the genome and genomic coverage of 13.9 GB, the average marker density was one marker per 0.95 Mb. The highest marker density with one marker every 0.54 Mb was observed on chromosome 2B, and the lowest density with one marker per 6.854 Mb on chromosome 4D (Fig. S1). The SNPs per chromosome ranged from 74 (4D) to 1486 (2B). The minimum and maximum SNPs detected on individual chromosomes ranged from 574 (4A) to 1293 (7A), 378 (4B) to 1486 (2B), and 74 (4D) to 255 (1D) at A, B, and D genome. The homoeologous group 7 had a maximum number of SNPs, 2838, followed by group 2 with 2652 SNPs, and group 4 had a minimum number of 1026 SNPs.

Genetic analysis of SNPs

The heterozygosity among all SNPs ranged from 0.00 to 0.149. Twenty-one SNPs showed high heterozygosity (>0.10), and 3,095 SNPs showed zero heterozygosity (Fig. 1B, Table S2). The heterozygosity was less than 0.05 for around 97% of the SNPs. SNPs in chr7A showed the highest heterozygosity, with SNP S7A_717968190 having heterozygosity of 0.149. The three ABLs, namely GID7396143, GID7399653, and GID7400318, were highly heterologous with heterozygosity equal to 0.16, 0.15, and 0.12, respectively (Table S3). The MAF ranged from 0.05 to 0.50 in the present study. Fifty-five SNPs were found with MAF =0.50 (where both the alleles of these SNPs were equally distributed across the panel). The average observed MAF was 0.23, while most ranged from 0.05 to 0.10, i.e., 2923 SNPs (Fig. 1C). The PIC values ranged from 0.085 to 0.250, with 54 SNPs having a PIC value equal to 0.250 (Fig. 1D). 5353 SNPs (36%) showed a range of 0.24−0.25 for PIC and were the highest among the SNPs.

Population structure analysis

The STRUCTURE results for the 141 ABLs showed that ΔK was highest at K = 2, indicating the presence of two major subgroups (Fig. 2A). The output of the STRUCTURE for K = 2 is given in Fig. 2B. The structure out showed very high ΔK for K = 2, compared to other K values, only a few genotypes represented complete distinction between the subgroups. Some genotypes across the subgroups showed near equal distribution of alternative alleles. Similar results were revealed by the PCA analysis where the genotypes of the subgroups could be seen clustered near each other when the first three principal components were observed (Fig. 2C). The dendrogram based on PCA also showed similar results (Fig. 3).

Analysis of molecular variance (AMOVA)

The analysis was performed considering K = 2. The variance among groups for K = 2 explained only 1% of the total variance and 99% variance within the subgroups. It indicates that the genetic differentiation within subgroups was high. Simultaneously, it was low among the subgroups (Table 2, Fig. 4), which may be due to the low genetic differentiation or excessive gene traffic between the population set. The high gene traffic is supported by the variation trends of both percentages of variance explained and haploid Nm values. The Nm values (haploid number of migrants) was very high, with Nm = 86.428.

Allelic pattern covering the populations

The mean Shannon’s information index (I), the diversity index (h), and the unbiased diversity index (uh) were high (0.477, 0.313, and 0.320), which indicated its high diversity (Table 3, Fig. 4). The mean of different alleles (Na) and the effective alleles (Ne) of the K = 2 population set were 1.977 and 1.519. The different alleles (Na) in the subgroup 1 (G1) and subgroup 2 (G2) were 1.954 and 2.0, respectively, indicating comparatively better diversity in G2 than G1 (Table 3).

Similarly, the numbers of effective alleles (Ne) was relatively higher in G2 (1.525) than in G1 (1.512). However, only G2 showed a significant percentage of private alleles, i.e., 4.6%. G1 contained 22 members with private allele ranging from 7 in GID7399640 to 149 in GID7400491, while G2 contained 90 members with private alleles ranging from 3 in GID7396143 to 188 in GID7400293 (Table 3 & Table S4). The diversity indexes I, h and uh for G1 & G2 were 0.468 and 0.487, 0.307 and 0.318, 0.320 and 0.321, respectively, indicating relatively higher diversity within G2 than G1.

Conclusion

In the present study, we applied GBS-based SNP to learn GBS-SNP markers’ usefulness for diversity analysis in 141 elite wheat breeding lines. Despite very designated, our advanced breeding lines panel was found to be genetically diverse, which could be instrumental for future South Asian breeding programs to develop new elite wheat varieties of alluring traits, i.e., high yield, biotic and abiotic resistance. Besides, the present study dappled two subgroups that were enlightened by their parentage and selection history. The low heterozygosity detected among elite advanced wheat breeding lines within subgroups and the moderate divergence among subgroups suggested that the elite advanced wheat breeding lines could be used further for GWAS studies.