Genome-wide association analysis unveils novel QTLs for seminal root system architecture traits in Ethiopian durum wheat

Abstract

Background

Genetic improvement of root system architecture is essential to improve water and nutrient use efficiency of crops or to boost their productivity under stress or non-optimal soil conditions. One hundred ninety-two Ethiopian durum wheat accessions comprising 167 historical landraces and 25 modern cultivars were assembled for GWAS analysis to identify QTLs for root system architecture (RSA) traits and genotyped with a high-density 90 K wheat SNP array by Illumina.

Results

Using a non-roll, paper-based root phenotyping platform, a total of 2880 seedlings and 14,947 seminal roots were measured at the three-leaf stage to collect data for total root length (TRL), total root number (TRN), root growth angle (RGA), average root length (ARL), bulk root dry weight (RDW), individual root dry weight (IRW), bulk shoot dry weight (SDW), presence of six seminal roots per seedling (RT6) and root shoot ratio (RSR). Analysis of variance revealed highly significant differences between accessions for all RSA traits. Four major (− log₁₀P ≥ 4) and 34 nominal (− log₁₀P ≥ 3) QTLs were identified and grouped in 16 RSA QTL clusters across chromosomes. A higher number of significant RSA QTL were identified on chromosome 4B particularly for root vigor traits (root length, number and/or weight).

Conclusions

After projecting the identified QTLs on to a high-density tetraploid consensus map along with previously reported RSA QTL in both durum and bread wheat, fourteen nominal QTLs were found to be novel and could potentially be used to tailor RSA in elite lines. The major RGA QTLs on chromosome 6AL detected in the current study and reported in previous studies is a good candidate for cloning the causative underlining sequence and identifying the beneficial haplotypes able to positively affect yield under water- or nutrient-limited conditions.

Background

Ethiopian farmers have grown tetraploid wheat (Triticum turgidum ssp. durum) since its introduction in the northern highlands of the country around 3000 BC [1]. Cultivation was mostly under adverse environmental conditions that likely favored the development of a broad gene pool of durum wheat landraces adapted to various environmental conditions. Ethiopian durum wheat landraces provide a rich and yet untapped native biodiversity [2]. Vavilov [3] and Zohary [4] reported the presence of high-genetic diversity in cultivated tetraploid wheat and recent studies highlighted the uniqueness of Ethiopian durum landraces from the Fertile Crescent collections (primary center of domestication) and considered Ethiopia as a possible second domestication center for the crop [5]. Previous studies, carried out with phenotypic [2, 6,7,8] and molecular approaches [9,10,11,12], have indicated Ethiopian durum germplasm to be a highly diverse and potentially unique source of valuable traits [13,14,15]. This is basically due to the wide range of agro-ecological conditions (altitude in a range of 1600 to 3000 masl) coupled with diverse farmers’ culture [9]. Notably, more than 7000 Ethiopian durum wheat landrace accessions are conserved in the Ethiopian Biodiversity Institute (EBI) gene bank [16]. However in recent time, durum wheat cultivation has been largely replaced by bread wheat varieties developed from international and national breeding programs throughout the country [17].

Roots play a key role in nutrient and water uptake, soil anchoring and mechanical support, storage functions, and as the major interface between the plant and various biotic and abiotic factors in the soil environment. Root system architecture (RSA) describes the shape and structure of the root system, both of which have great functional importance [18, 19] and plays a pivotal role in crop performance, especially for cultivation under non-optimal nutritional and water source conditions [20,21,22]. Due to recurrent climate change, declining of soil fertility and water availability, enhancing the genetic capacity to capture the available soil resources is considered a primary target for breeding resource-use efficient crops [20, 23, 24]. Hence, RSA has been an active research topic for the last couple of decades and since then different RSA ideotypes have been proposed and investigated in crops [25,26,27]. The narrow-and-deep or wide-and-shallow root ideotypes have been studied for their effects in nutrient acquisition and drought resistance in crops [28,29,30,31]. Deep and narrow-angled roots could allow plants to exploit more effectively water and nitrogen that are often found in deeper soil layers [29, 30, 32], while shallow wider angled roots enable plants to more effectively uptake nutrients such as phosphorous that are abundantly found at shallower depths in the soil [33].

The genetic basis of RSA traits in durum wheat has been investigated with both linkage and association mapping using durum wheat recombinant inbred line (RIL) populations and/or elite durum wheat panels suitable for association mapping [19, 21, 34,35,36,37]. This notwithstanding, beside the recent studies by Roselló et al. [38] and Ruiz et al. [39], durum wheat landraces have not been extensively studied so far. Ethiopian durum wheat landraces are particularly rich in genetic diversity and thus are very valuable to dissect the genetic basis of governing the variability of RSA traits. Hence, this study aimed to conduct a genome-wide association analysis for root system architecture traits in Ethiopian durum wheat comprising historical landraces (167) and modern cultivars (25) to identify RSA quantitative trait loci (QTLs) of potential interest for marker-assisted selection.

Results

Phenotypic variation among RSA traits

A total of 2880 seedlings and 14,947 seminal roots were processed and measured for various RSA traits (Additional file 2: Table S2). Analysis of variance (ANOVA) for the studied RSA traits is presented in Table 1.

Table 1 ANOVA and heritability results for the root system architecture traits measured in 12-day-old seedlings of 192 Ethiopian durum wheat accessions

The ANOVA results indicate the presence of highly significant variation among accessions for all RSA traits. In particular, the seminal root angle ranged from 45.7 to 130.5° with a mean value of 97.3° while the total and average root length and number of roots ranged from minimum values of 66.2 cm, 16.5 cm and 3.4 to maximum values of 195.4 cm, 36.9 cm and 6.7, respectively. The root and shoot dry weight varied from minimum values of 27.7 and 34.7 g to maximum values of 115.0 and 116.6 g, respectively. The coefficient of variance (CV) of RSA traits ranged from 8.38% for average root length (ARL) to 14.63 for root growth angle (RGA). Individual root dry weight (IRW) and bulk root dry weight (RDW) also scored high CV, with a value of 14.55 and 14.22%, respectively. The frequency distribution of most RSA traits was normal except for RT6 that showed a bi-modal distribution (Fig. 1).

Fig. 1

Distribution frequency for RSA traits measured from 12-day-old seedlings in 192 Ethiopian durum wheat accessions. See Table 6 for trait abbreviations

Most RSA traits showed high level of broad sense heritability (H²). Bulk root dry weight (RDW), average root length (ARL) and bulk shoot dry weight (SDW) showed the top three values (91.3, 91.0 and 90.4%, respectively) while the presence of the 6th root showed the lowest value (67.0%).

Correlation among RSA traits

Several strong correlations were observed between RSA traits (Fig. 2). Highly significant positive correlations were detected for RDW vs. IRW (0.93), RDW vs. SDW (0.92) and IRW vs. SDW (0.84). Strong correlations were recorded between TRN and RT6; TRL and ARL with a correlation coefficient of 0.84 and 0.82, respectively. The initial thousand grain weight showed no significant correlation with any RSA trait suggesting that variation of RSA traits did not have maternal etiology caused by variation in seed size.

Fig. 2

Correlation coefficient and level of significant for the initial thousand grain weight and RSA traits measured in 12-day-old seedlings of 192 Ethiopian durum wheat accessions

Landraces showed a wider range of variability than cultivars in most RSA traits although the latter outperformed the former for some traits (Table 2 and Additional file 9: Figure S2.). For instance, the cultivars mean values for root and shoot dry weight were 90.3 and 92.8 mg, while landraces scored only 56.9 and 66.5 mg for the same traits, respectively. Cultivars also performed better than landraces for TRL and ARL while TRN and RT6 were the only two RSA traits for which landraces showed slightly higher mean values than cultivars.

Table 2 Mean and range values of 25 cultivars and 167 landraces for RSA traits

Population structure and linkage disequilibrium decay analysis

According to population structure analysis, the panel was subdivided into three subpopulations of 75, 27 and 90 accessions each (Fig. 3a, b and Additional file 3: Table S3). All 26 cultivars clustered into subpopulation 2 except for ‘Selam’ that grouped in subpopulation 1. Clustering analysis indicated that SNP data failed to group landraces clearly based on their geographical backgrounds and accessions were admixed into the three subpopulations irrespective of their geographic origin. Box plot of the three sub-populations inferred from STRUCTURE analysis for the mean values of RSA traits is reported in Additional file 9: Figure S3.

Fig. 3

Population structure and kinship-matrix similarity analysis for 192 Ethiopian durum wheat accessions. Heat-map clustering results based on the kinship matrix from tag-SNP (r² = 1) by identity-by-state (IBS) algorithm (a). Population structure plot and K1, K2 and K3 represents subpopulations 1, 2 and 3, respectively (b). The black-dash lines separated the panel into three subpopulations. Accessions arrangement was based on the order of heat-map kinship result. The color represents the membership of each accession in the STRUCTURE-inferred subpopulations. The color of the legend indicates the level of kinship similarity of the heat-map

The mean genome wide r² value was 0.12, with 55% of the pair-wise linkage disequilibrium comparisons showing significant association at P < 0.01. Chromosome 3B scored the highest mean value (r² = 0.19) with 64% significant pair-wise LD comparisons. On the other hand, 7A scored the lowest mean r² value (0.11) and 48% of pairwise LD comparisons were significant. The genome-wide LD decayed below r² = 0.3 (the standard critical threshold) at 2.25 cM. This defines the ±2.25 cM as the genome-wide critical distance to detect linkage and, therefore, as the QTL confidence interval around the QTL-tag SNP, i.e. the SNP found at the peak of the corresponding QTL. The specific critical r² value beyond which LD is due to true physical linkage was 0.15 and the intersect of the threshold with the LD decay curve was at 5.75 cM.

GWAS analysis of RSA traits

After filtering SNP data and following imputation, a total of 10,789 polymorphic SNP markers (4591 and 6198 SNPs from A and B genomes, respectively) were used for marker-traits association (MTA) analyses. The mixed linear model with population structure and kinship matrix was chosen for MTA analysis, as the quantile-quantile (Q-Q) plot showed that the observed MTA P-values were close to the expected distribution (Additional file 9: Figure S4). A total of 275 QTLs with various significant values were identified for the tested RSA traits. The only four major QTLs above the experiment-wise threshold (− log₁₀P ≥ 4) were EPdwRGA-6A, EPdwRDW-4A, EPdwiTGW-3B.1 and EPdwIRW-5A with values of 6.85, 4.34, 4.15 and 4.06 which accounted for 16.08, 8.41, 8.71 and 8.03% of the phenotypic variation, respectively. Thirty-four QTLs reached the marker-wise threshold of – log₁₀P ≥ 3 in which the highest number was identified for TRN with eight QTLs followed by SDW and IRW each with six nominal QTLs. Additionally, three nominal QTLs were identified for TRL, iTGW and RT6, two for RDW and only one for RGA, ARL and RSR. The other 237 QTLs with a marker-wise threshold of – log₁₀P ≥ 2 were identified as suggestive QTLs. The major and nominal QTLs are reported in Table 3 while the complete list of identified QTLs with the marker-wise threshold value of –log₁₀P ≥ 2 are reported in Additional file 4: Table S4. Thirteen markers showed significant associations for more than one RSA trait that could be due to either a pleiotropic effect or tight linkage, hence considered as separate QTLs for corresponding traits (Table 4). Notably, the root growth angle QTL showed limited overlap with QTLs of other RSA traits.

Table 3 List of major and nominal QTLs for RSA traits identified in 192 Ethiopian durum wheat accessions

Table 4 Markers with a significant association/concurrent effect on more than one RSA trait

QTL clusters for RSA traits

The identified QTLs were further grouped into 15 RSA QTL clusters plus one distinct RGA QTL cluster on chromosome 6AL. Clustering was based on the significance of each QTL and its effects on various traits in this study and overlapping with QTLs from previously reported studies in bread and/or durum wheat (Table 5). Based on these criteria, a total of 103 QTLs were included in 16 QTL clusters. Cluster pairs were identified on chromosomes 1A, 3B and 7A while chromosomes 1B, 2A, 2B, 3A, 4A, 4B, 5A, 5B, 6A and 6B each harbored a single QTL cluster (Fig. 4a, b and Additional file 9: Figure S5).

Table 5 Main RSA QTL clusters identified in 192 Ethiopian durum wheat accessions and other studies

Fig. 4

Genetic map of RSA QTLs identified in 192 Ethiopian durum wheat accessions along with previously published studies projected onto SNP-based tetraploid consensus map published in Maccaferri et al. (2015). RSA QTL identified in the present study are listed on the left side of the chromosomes with their significance level: ** = marker-wise significance of P ≤ 0.01 (− log₁₀P ≥ 2); *** = marker-wise significance of P ≤ 0.001(− log₁₀P ≥ 3); **** = experiment-wise significance of P ≤ 0.05/marker-wise significance of P ≤ 0.0001 (− log₁₀P ≥ 4). RSA QTLs identified in previous studies, orange-filled bars for durum wheat and blue-filled bars for bread wheat, listed on the right side and references given in parentheses. Grey-filled bands are for RSA QTL clusters on chromosomes 1A and 1B (a) and a distinct root growth angle (RGA) QTL cluster identified on chromosome 6A from 105 to125 cM (b). Black-filled bars are for QTLs with R² < 5%; red bars for R² values from 5 to 10% and yellow bars for R² > 10%. The length of bars indicates the confidence interval of each QTL or QTL cluster. Manhattan plot for the major RGA QTL identified on chromosome 6AL (c)

QTL for seminal root length and number

EPdwTRL-1B, EPdwTRL-4B and EPdwTRL-5A were the three nominal QTLs identified for TRL on chromosomes 1B (Fig. 4a), 4B and 5A, respectively. Other suggestive TRL QTLs were identified on all chromosomes except for chromosome 6A. For ARL, only one nominal QTL (EPdwARL-2A) was detected on chromosome 2A, while other suggestive QTLs were detected for across all chromosomes. Seven nominal QTLs were detected for TRN: three (EPdwTRN-4B.1, EPdwTRN-4B.2 and EPdwTRN-4B.3) were mapped on chromosome 4B, two (EPdwTRN-1A.1 and EPdwTRN-1A.2) on chromosome 1A (Fig. 4a) and the other two (EPdwTRN-1B and EPdwTRN-7A) on chromosomes 1B and 7A, respectively. For the presence of the sixth seminal root, three nominal QTLs (EPdwRT6-4B.1, EPdwRT6-4B.2 and EPdwRT6-4B.3) were mapped on chromosome 4B (Table 3). The allelic distribution and frequency of TRN and TRL QTL-tagging SNPs with phenotypic effect (R²) > 5% are reported in Additional file 6: Table S6 and Additional file 7: Table S7, respectively.

QTL for seminal root growth angle

The QTL with the largest effect (R² = 0.16) on RGA (EPdwRGA-6A) was identified on chromosome 6A. Within the confidence interval of this QTL, six SNPs (IWB35245, IWB71122, IWB24306, IWB57413, IWB10077 and IWB74235) showed significant effects for the trait (Fig. 4c; Additional file 4: Table S4). The confidence interval of this major RGA QTL (from 105 to 125 cM) overlapped with the confidence interval of RSA QTLs previously reported in the same region (Fig. 4b). Other suggestive RGA QTLs were identified on chromosomes 1A, 2B, 3A, 3B, 4A, 5B, 6B, 7A and 7B (Additional file 4: Table S4). Notably, RGA QTLs showed no clustering with other RSA QTLs. The allelic distribution and frequency of RGA QTL-tagging SNPs with phenotypic effect > 5% is reported in Additional file 5: Table S5.

QTL for root and shoot dry weight

Two major QTLs (EPdwRDW-4A and EPdwIRW-5A) were identified for bulk and individual root dry weight on chromosomes 4A and 5A, respectively. Two nominal QTLs were identified for RDW (EPdwRDW-1B and EPdwRDW-3A) on chromosomes 1B and 3A. As to individual root weight six nominal QTLs (EPdwIRW-1B, EPdwIRW-2B, EPdwIRW-5B.1, EPdwIRW-5B.2, EPdwIRW-6B and EPdwIRW-7A) were identified on chromosomes 1B, 2B, 5B (two QTLs), 6B and 7A, respectively. Six nominal QTLs (EPdwSDW-1A, EPdwSDW-1B, EPdwSDW-3B, EPdwSDW-4A, EPdwSDW-4B and EPdwSDW-5B) were identified for SDW. The QTLs for these three traits repeatedly clustered nearby or in single QTLs (Table 5). The allelic distribution and frequency of IRW QTL-tagging SNPs with phenotypic effect > 5% is reported in Additional file 8: Table S8.

Discussion

In the present study, 12-day-old seedlings of 192 Ethiopian durum wheat accessions, predominantly landraces, were phenotyped in controlled conditions to identify the root system architecture (RSA) QTL through GWAS analysis. Moderate to high heritability values, ranging from 67 to 91%, were recorded for all RSA traits, confirming them as potential targets for wheat improvement.

The linkage disequilibrium analyzed from 10,789 polymorphic SNPs indicated that LD decays to the threshold value of r² = 0.3 (the generally accepted limit to detect association with a QTL) at 2.25 cM that was in agreement with the LD decay value previously detected by Liu et al. [15]. Maccaferri et al. [40, 41] specified the LD decays at 2.20 cM for the panel comprising 183 elite durum wheat cultivars and lines from Mediterranean countries, the Southwestern USA and Mexico.

The RSA QTL-clusters included either single loci with concurrent effects on different RSA traits or tightly linked loci not resolved by recombination [42], most of which overlapped with previously identified RSA QTL clusters. QTL mapping for RSA traits of wheat based on designed bi-parental populations was recently reviewed by Soriano and Alvaro [43] compiling the results of 27 bread and three durum wheat studies for a total of 754 QTLs.

Root length and number at the seedling stage are potential candidates for marker-assisted breeding applications aimed at enhancing early rooting capacity [21]. One novel QTL for TRN, EPdwTRN-4A, was discovered in the present study on the short arm of chromosome 4A. The other TRN QTL identified on the short arm of chromosome 1A overlaps with the TRN QTL reported by Maccaferri et al. [21]. The confidence interval of the TRN QTL on the short arm of chromosome 1B overlapped with the confidence interval of the TRN QTL identified by Christopher et al. [44] and under the 8th root metaQTL (Root_MQTL_8) reported by Soriano and Alvaro [43]. Other nominal TRN QTL identified on the short arm of chromosome 4B overlapped with TRN QTL reported by Ren et al. [45]. The other two TRN QTLs detected on the long arm of chromosome 4B and short arm of chromosome 7A both overlapped with a TRN QTL reported in Maccaferri et al. [21]. Chromosome 4B showed three strong QTLs (EPdwRT6-4B.1, EPdwRT6-4B.2 and EPdwRT6-4B.3) for the development of more than five seminal roots per plantlet.

For root length, the other important trait, d three nominal QTLs were identified for TRL and one for ARL. One novel QTL for TRL, EPdwTRL-4B, was mapped on the long arm of chromosome 4B. The TRL QTL identified on the short arm of chromosome 1B overlaps with TRL QTL reported by Petrarulo et al. [36] and Liu et al. [46] and the other one detected on the telomeric region of chromosome 5A overlapped with a TRL QTL reported by Maccaferri et al. [21]. The nominal ARL QTL (EPdwARL-2A) identified on chromosome 2A with a concurrent effect on TRL, SDW, RDW and IRW, is novel since it was not reported in any of the previous studies considered for this meta-analysis based on the tetraploid consensus map.

Among the other essential RSA traits, as to root growth angle (RGA), a pivotal trait influencing RSA and its functions, the most notable QTL (EPdwRGA-6A) was identified on the long arm of chromosome 6A, similarly tothat reported by Maccaferri et al. [21], QRga.ubo-6A.2, using 183 elite cultivars and lines representing the main breeding pools from Mediterranean countries (particularly ICARDA and Italy), the Southwestern USA and CIMMYT. Additionally, Alahmad et al. [47] recently reported sizeable and highly significant effects on RGA of the same region of chromosome 6AL. The concomitant effects of the chromosome 6AL on RGA observed in widely different germplasm pool underline the importance of further studies to better characterize the effects of the different haplotypes present at this major QTL. Notably, a novel nominal RGA QTL (EPdwRGA-4A) was detected on the long arm of chromosome 4A.

An additional novel major RDW QTL (EPdwRDW-4A) with concurrent effects on SDW, TRN and TR6 was mapped on the short arm of chromosome 4A. A novel RDW QTL (EPdwRDW-3A) was also identified on the long arm of chromosome 3A. EPdwSDW-3B and EPdwSDW-4A were the two newly discovered nominal SDW QTLs on the short arm of chromosome 3B and long arm of chromosome 4A, respectively. Four novel IRW QTLs (EPdwIRW-5B.1, EPdwIRW-5B.2, EPdwIRW-6B, EPdwIRW-7A) were discovered on the short arm of chromosome 5B (the first two), long arm of chromosome 6B and short arm of chromosome 7A, respectively. Iannucci et al. [37] noted the absence of a clear relationship between plant height and root development and added diverse and controversial speculations from a number of previous studies which are probably due to the different conditions and growth stages in which the root traits were evaluated. Some authors reported different genetic control between shoot and root growth [35, 48, 49] while others have reported a negative correlation [50]. Bai et al. [51] investigated a set of NILs for a number of Rht loci/alleles and showed clear effects on both shoot and root traits.

Conclusions

Among the four major and 34 nominal RSA QTLs identified in the current study, 14 are novel, hence showing the suitability of Ethiopian landraces for studies aimed at the dissection of the QTL and the identification of novel haplotypes. The remaining 20 RSA QTLs concomitantly identified in this and previous studies provide valuable information on their role across diverse genepools, an important prerequisite to prioritize QTLs for marker-assisted selection aimed at enhancing crop productivity based on the use of RSA traits as proxies. A cluster of RGA QTLs was identified on the long arm of chromosome 6A with a major QTL (EPdwRGA-6A) with a notable phenotypic effect on RGA (R² = 0.16). This result coupled with those reported in previous RSA studies [21, 47] highlights and reinforces EPdwRGA-6A as a strong candidate for further studies aimed at cloning the causative sequences and identifying the beneficial haplotypes able to positively affect yield under water- or nutrient-limited conditions.

Methods

Plant materials

One hundred ninety-two Ethiopian durum wheat accessions were used to assemble the GWA mapping panel. The collection included 167 landraces and 25 cultivars collected and maintained as single seed descent (SSD) progenies at the Debre Zeit Agricultural Research Center (DZARC) and Sinana Agricultural Research Center (SARC) in Ethiopia.

Landrace collections were originally collected from major wheat-producing areas of Ethiopia, including Bale, Gondar, Gojjam, Shewa, Tigray and Wollo. Twelve Ethiopian durum wheat landraces currently cultivated in the USA are included in the panel. Cultivars were released in the years between 1994 and 2010 from DZARC and SARC and have been/are being cultivated in Ethiopia. Details of accessions used for the current study are summarized in Additional file 1: Table S1.

Root system architecture phenotyping

Seminal RSA traits were characterized using the protocol described by Canè et al. [19] and later used by Maccaferri et al. [21] with minor adjustments in the present work. Seeds were first weighed to measure thousand grain weight that was later used as a covariate in order to account for maternal effects on RSA traits due to seed size. Twenty seeds per accession were treated in 0.15% Panoctine solution and dried before pre-germinating them in Petri dishes on wet-filter-paper at 28 °C for 24 h. Then, five similar seeds with homogenous seminal root emission were positioned 7-cm apart on a wet-filter-paper sheet moistened with distilled water and placed on a vertical black rectangular (42.5 × 38.5 cm) polycarbonate plate for root obscuration.

Root traits were then measured in plantlets grown in a growth chamber for 12 days at 22 °C (day)/18 °C (night) under a 16-h photoperiod and light intensity of 400 μmol m^− 2 s^− 1 photosynthetically active radiation (PAR). The experiment was conducted adopting a randomized complete block design (RCBD) with three independent replications grown in the growth chamber. The experimental unit included five homogenous seedlings of each accession and hence one screening plate corresponded to one genotype. Blocking was introduced to control for possible differences in growth rate and normalization of the blocking effect (linear adjustment, whenever significant) was undertaken. Due to the high number of genotypes under evaluation and the time required for root preparation and root image acquisition, genotypes were divided into sets of 25–30 accessions that were considered as blocks. Blocks included accessions phenotyped at the same date and kept on shelves in the growth chamber that are positioned at the same distance from the floor under uniform light conditions (see Additional file 9: Figure S1).

Data for the following RSA traits were taken based on single-plantlet basis (Table 6): root growth angle (RGA) measured as the linear distance between the two most external seminal roots of each plantlet at 3.5 cm from the seed tip and then converted to degrees (Fig. 5a, b); total root length (TRL); average root length (ARL); total root number (TRN); presence of six seminal roots (RT6). Total root length and root growth angle were measured on plantlet images (Fig. 5c) using GIMP (GNU Image Manipulation Program) and ImageJ [52]. Average root length was estimated as total root length divided by total root number. Bulked roots and shoots from each experiment were cut and dried in an oven for 48 h to measure root dry weight (RDW) and shoot dry weight (SDW), respectively. Individual root dry weight (IRW) was derived from the result of the bulk root dry weight divided by the total root number that could be used as a proxy to measure root thickness.

Table 6 Summary of acronyms used for root system architecture (RSA) traits and their measuring unit

Fig. 5

Root growth angle of seminal roots in 12-day-old seedlings of ‘Gondar’ landrace with narrow growth angle (a) and ‘Obsa’ cultivar with wide growth angle (b) measured as the linear distance (red segment) of the two most external roots (green segments) at 3.5 cm from the tip of the seed and later converted into degrees. Example of a root sample ready for image capturing for further root length and root growth angle measurement (c)

Phenotypic data analysis

Analysis of variance (ANOVA) was conducted including replications, blocks and accessions. Block effect was controlled using the mean of each set of genotypes included in the same block and used to correct the corresponding single values, whenever significant, with a linear regression method. The weight of each individual seed was used as a covariate to correct for any possible variation caused by maternal effects. In addition, the trait was subjected to GWA analysis along with other RSA traits.

Broad sense heritability (H²) of RSA traits was calculated with the mean values of each experiment among the three replications according to the formula:

$$ {H}^2=\frac{\sigma^2g}{\sigma^2g+{\sigma}^2e/r} $$

Where σ²_g (genetic variance) was calculated as (MS_genotypes – MS_residual)/r; σ²_e (the residual variance) = MS_residual, r the number of replications and MS the mean square value. The coefficient of variance (CV) was calculated for all RSA traits except for the presence of the 6th root, the only trait with discrete values.

Genotypic data and imputation

A pooled tissue sample of 25 one-week-old plantlets, from the same seed source used to phenotype RSA traits, was used for genomic DNA extraction for each accession. DNeasy 96 Plant Kit (Qiagen GmbH, Hilden, Germany) was used to extract the genomic DNA. Genotyping was done with the high-density Infinium® iSelect® Illumina 90 K wheat SNP array [53] and SNP calling and clustering were made with the GenomeStudio v2011.1 software (Illumina, San Diego, CA, USA). Calls showing residual heterozygosity were assigned as a missing value. SNP markers with < 0.05 minor allele frequencies (MAF) and markers with > 0.1 missing values per accession were excluded. After filtering, imputation of the missing data was computed using Beagle 4.0 [54]. Owing to the high level of homozygosity, imputation disregarded any phased reference populations. Twenty-five markers were considered in the imputation rolling window (twice the average number of marker present in a 5 cM interval), with an overlap of a single marker, the typical number of markers included in a 0.5 cM interval. Since imputation accuracy was not improved by using other parameters, default values were kept.

The high-density consensus map of tetraploid wheat generated by Maccaferri et al. [41] was used to identify chromosome positions of SNPs and markers with unknown positions were removed.

Population structure and kinship analysis

For population structure analysis, a Bayesian model-based (Markov Chain Monte Carlo) clustering approach was used in STRUCTURE v.2.3 [55]. Haploview v4.2 [56] “Tagger” function (based on analysis of marker pairwise r² values) was used to select tag-SNPs for population structure analysis with a tagger filter set at r² = 0.5 and 1496 tag-SNPs were selected.

To infer the optimal sub-populations number, an ad hoc quantity (∆K) was calculated based on the second order rate of change of the likelihood (Evanno et al., 2005) and in this analysis approach, the ∆K shows a clear peak at the ideal number of sub-populations. To perform this, 10 sub-populations with 20 independent iterations for each sub-population were done adopting an admixture model of population structure with correlated allele frequencies and 50,000 lengths burn-in period and 100,000 Markov Chain Monte Carlo (MCMC) replications after burn-in were applied for each iteration.

Additionally, the Haploview “Tagger” function was used to select tag-SNPs for kinship matrix (K) analysis with a tagger filter set at r² = 1 and 4842 tag-SNPs were selected, calculated in TASSEL v.5.2 [57] and incorporated in the mixed linear model (MLM) along with the population structure (Q) value for GWAS analysis.

Linkage disequilibrium (LD) and GWAS analysis

The LD r² values between pairwise intra-chromosomal SNPs were calculated with TASSEL v.5.2 and LD decay curve was fitted by a smoothing spline regression line at the genome level according to Hill and Weir function [58] in r environment [59]. The specific critical r² value beyond which LD is due to true physical linkage was determined by taking the 95th percentile of r² data of unlinked marker pairs [60]. In order to control the rate of false-positive associations, a MLM model [61] with population structure and kinship covariates was applied for the GWAS analyses. Hence, all SNP markers and the phenotypic data generated for the nine RSA traits were used to conduct the MTA analysis.

Three levels of significance were introduced according to Maccaferri et al. [21] for reporting the GWAS-QTLs: (i) experiment-wise P ≤ 0.05 (marker-wise P ≤ 0.0001, − log₁₀P ≥ 4) for “major QTLs”; (ii) marker-wise P ≤ 0.001 (− log₁₀P ≥ 3) for “nominal QTLs”; (iii) marker-wise P ≤ 0.01, (− log₁₀P ≥ 2) for “suggestive QTLs”. The experiment-wise threshold was established according to the number of ‘independent SNP tests’ that was estimated in Haploview using the tagger function of r² = 0.3 [62] and the total number (816) of tag-SNPs. Bonferroni test adjusted for multiple marker tests (P ≤ 0.05) was equal to – log₁₀P = 4.21 (rounded to 4.00). Hence the experiment-wise, Bonferroni-corrected significance threshold at P = 0.05 matched to a marker-wise threshold of – log₁₀P ≥ 4. Significance intervals of identified QTLs were reported as the intervals after including all SNPs associated with the trait with P ≤ 0.01 (marker-wise) and in LD of r² ≥ 0.3. Confidence intervals were defined based on the GWAS-QTL peak ±2.25 cM on both map sides.

The relative positions of RSA QTLs identified in this study along with other previous studies [14, 21, 34, 36, 37, 44,45,46, 51, 63,64,65,66,67] were compared based on the projected QTL peaks and confidence intervals on the tetraploid wheat consensus map [41].