Exploring the Genetic Architecture of Root-Related Traits in Mediterranean Bread Wheat Landraces by Genome-Wide Association Analysis

Abstract

Background: Roots are essential for drought adaptation because of their involvement in water and nutrient uptake. As the study of the root system architecture (RSA) is costly and time-consuming, it is not generally considered in breeding programs. Thus, the identification of molecular markers linked to RSA traits is of special interest to the breeding community. The reported correlation between the RSA of seedlings and adult plants simplifies its assessment. Methods: In this study, a panel of 170 bread wheat landraces from 24 Mediterranean countries was used to identify molecular markers associated with the seminal RSA and related traits: seminal root angle, total root number, root dry weight, seed weight and shoot length, and grain yield (GY). Results: A genome-wide association study identified 135 marker-trait associations explaining 6% to 15% of the phenotypic variances for root related traits and 112 for GY. Fifteen QTL hotspots were identified as the most important for controlling root trait variation and were shown to include 31 candidate genes related to RSA traits, seed size, root development, and abiotic stress tolerance (mainly drought). Co-location for root related traits and GY was found in 17 genome regions. In addition, only four out of the fifteen QTL hotspots were reported previously. Conclusions: The variability found in the Mediterranean wheat landraces is a valuable source of root traits to introgress into adapted phenotypes through marker-assisted breeding. The study reveals new loci affecting root development in wheat.

1. Introduction

Wheat is the most widely cultivated crop in the world, covering around 219 million ha (Faostat 2017, http://www.fao.org/faostat/). It is a staple food for humans, as it provides 18% of daily human intake of calories and 20% of protein (http://www.fao.org/faostat/). Global wheat demand is estimated to increase by 60% by the year 2050 [1], so wheat production will need to rise by 1.7% per year until then. Achieving this objective is a great challenge under the current climate change scenario, as the prediction models estimate a precipitation decrease of 25% to 30% and a temperature increase of 4 °C to 5 °C for the Mediterranean region [2]. It is well known that wheat production is greatly affected by environmental stresses such as drought and heat [3] that negatively affect yield and grain quality [4]. Drought is considered the greatest environmental constraint to yield and yield stability in rainfed production systems [5]. Environmental effects on yield in the Mediterranean Basin have been estimated at 60% for bread wheat [6] and 98% for durum wheat [7]. The expected effects of climate change and the declining availability of water and chemical fertilizers will require the release of cultivars with an enhanced genetic capacity to maintain acceptable yield levels and yield stability under harmful environmental conditions [8,9]. To cope with the challenges of climate change, breeders are particularly challenged to stretch the adaptability and performance stability of new cultivars, so many improvement programs are focussing on breeding for adaptation [10].

Plants respond and adapt to water deficit using various strategies that have evolved at several levels of function and are components of the conceptual framework developed by Reynolds et al. [11], which defines drought resistance in terms of dehydration escape, tolerance, and avoidance. Traits defining root system architecture (RSA) are critical for wheat adaptation to drought environments and non-optimal nutritional supply conditions [12]. Besides, water-use efficiency (WUE) can be significantly increased by optimizing the anatomy and growth features of roots [13]. Root traits are critical for drought tolerance due to its role in plant performance and the acquisition of nutrients and water from dry soils [14]. The wheat plant includes two types of roots: seminal (embryonal) and nodal (crown or adventitious or adult root system). The seminal roots are the first to penetrate the soil and remain functional during the whole plant cycle [9,15]. A correlation between seminal and adult roots in terms of size, dry-weight, or even specific architectural features have been reported [9,13]. Since the evaluation of RSA features in the field is very difficult, expensive, and time-consuming when a large number of genotypes need to be phenotyped, several studies have been carried out at early growth stages to allow an optimal screening of RSA traits [8,12,16,17,18]. Maccaferri et al. [9] observed that among RSA traits, those involving the root structure and related to the uptake of nutrients and water are root length, surface area and volume, and the number of roots, while root diameter is significantly associated with drought tolerance. Another RSA trait of interest in wheat is the seminal root angle (SRA), whose features suggest that narrow angles could lead to deeper root growth to obtain water from deeper soil layers and hence maintain higher yields [5,13].

Identifying quantitative trait loci (QTLs) and applying marker-assisted selection is of particular interest for RSA because the trait is important but difficult to phenotype. In the last few years, genome-wide association studies (GWAS) have become very popular because of their use of germplasm collections with wider variability than the classical bi-parental crosses. These collections allow many recombination events to be detected, making the association between genotype and phenotype more accurate. Collections of landraces are an ideal subject of GWAS [19] since they are genetically diverse repositories of unique traits that have evolved in local environments characterized by a wide range of biotic and abiotic conditions. Several studies have shown that Mediterranean wheat landraces possess a wide genetic background for root architecture, yield formation, stress tolerance, and quality traits [17,18,19,20,21,22]. In the current study, a GWAS for three RSA traits and two related traits was performed on a panel of 170 bread wheat (Triticum aestivum L.) landraces from 24 Mediterranean countries with the following goals: (1) to detect differences in RSA among genetic subpopulations previously distinguished in the panel, (2) to identify correlations among RSA and grain yield under rainfed conditions, and (3) to identify molecular markers and candidate genes linked to root-related traits and candidate gene models for the associations.

2. Materials and Methods

2.1. Plant Material

A germplasm collection of 170 bread wheat (Triticum aestivum L.) genotypes from the MED6WHEAT IRTA panel described by Rufo et al. [23] was used in this study. The panel was genotyped and characterized using the Illumina Infinium 15K Wheat SNP Chip at Trait Genetics GmbH (Gatersleben, Germany), and markers were ordered according to the SNP map developed by Wang et al. [24]. The collection was previously structured into three subpopulations (SPs) matching their geographical origin [23]: western (SP1, WM), northern (SP2, NM), and eastern Mediterranean (SP3, EM) (Supplementary Materials, Table S1). Additionally, the cultivars ‘Arthur Nick’, ‘Anza’, ‘Soissons’, and ‘Chinese Spring’ were included as checks.

2.2. Root Morphology and Statistical Analysis

Root analysis was performed following the protocol described by Canè et al. [8], which was slightly modified in the current study (Figure 1). Ten representative seeds were randomly chosen from each genotype, weighed, sterilized in a 10% sodium hypochlorite solution for 5–10 min, washed thoroughly in distilled water and placed on hydrated filter paper in a 140 mm Petri dish at 28 °C for 24 h. Subsequently, five seedlings were selected on the basis of a normal seminal root emergence and were spaced 8 cm from each other on a filter paper sheet placed on a vertical black rectangular polycarbonate plate (42.5 × 38.5 cm). Finally, each plate was covered with another wet sheet of filter paper. Distilled water was used for the plantlets’ growth. The plantlets were grown in a growth chamber for 14 days at 22 °C under a 16-h light photoperiod. In addition to the ten seed weight (SW), four other traits were scored for each genotype: total root number (TRN), shoot length (SL) from the seed to the tip of the longest leaf and SRA, obtained using a digital camera following the methodology described in Canè et al. [8]. The images were processed with ImageJ software [25]. The angle between the two external roots of each plantlet was measured at a distance of 3.5 cm from the tip of the seed. Finally, the roots were desiccated at 70 °C for 24 h to obtain the root dry weight (RDW).

The experimental design followed a randomized complete block with two replications in time. Means of five observational units for each genotype were used for TRN, RDW, and SL, while only three observational units were used for SRA because the two external ones were considered as border plantlets for root angle.

2.3. Grain Yield

Field experiments were carried out in 2016, 2017, and 2018 harvesting seasons in Gimenells, Lleida, north-east Spain (41°38’ N and 0°22’ E, 260 m a.s.l) under rainfed conditions. The experiments followed a non-replicated augmented design with two replicated checks (the cultivars ‘Anza’ and ‘Soissons’) and plots of 3.6 m². The experimental design is shown in Supplementary Materials, Figure S1. Sowing density was adjusted to 250 germinable seeds m². Weeds and diseases were controlled following standard practices at the site. The anthesis date was determined in each plot. Grain yield (GY, t ha⁻¹) was determined by mechanically harvesting the plots at maturity and expressed on a 12% moisture level.

2.4. Statistical Analysis

Phenotypic data for GY was fitted to a linear mixed model with the check cultivars as fixed effects and the row number, column number and cultivar as random effects following the SAS PROC MIXED procedure:

y = Xβ + Zγ + ε

(1)

where β is an unknown vector of fixed-effects parameters with known design matrix X, γ is an unknown vector of random-effects parameters with known design matrix Z, and ε is an unknown random error vector whose elements are no longer required to be independent and homogeneous.

Restricted maximum likelihood was used to estimate the variance components and to produce the best linear unbiased predictors (BLUPs) for the traits of each cultivar and year with the SAS-STAT statistical package (SAS Institute Inc, Cary, NC, USA).

Analyses of variance (ANOVA) were performed for the root traits, considering the genotypes and the replication as random effects in the model. Additionally, for a subset of 55 of the 141 structured landraces, selected as having an SP membership q > 0.8 (WM, 17; NM, 15; EM, 23), the sum of squares of the cultivar effect in the ANOVAs was partitioned into differences between SPs and differences within them. ANOVA for grain yield was performed for the complete collection, considering genotype, year, and the combination of genotype and year the sources of variation. Least squares means were calculated and compared using the Tukey HSD test at P < 0.05. Pearson correlation coefficients among root traits were computed. Repeatability (H) was calculated on a mean basis across two replications following the formula described by Harper [26] r = (B − W) / (B + ( (n −1 ) W)), where n is the number of genotypes and B and W the two variances from the ANOVA table: between (B) and within (W). Frequency distributions, ANOVAs, the Tukey test, and the Pearson correlation coefficients were calculated using the JMP v13.1.0 statistical package (SAS Institute Inc, Cary, NC, USA).

2.5. Genome-Wide Association Analysis

A GWAS was performed for the mean of measured root traits and from the BLUPs for GY per year and across years with TASSEL 5.0 software [27]. A mixed linear model (MLM) was conducted using the information of the genetic structure reported in Rufo et al. [23] as the fixed effect and a kinship (K) matrix, calculated using Haploview [28], as the random effect (Q + K model) at the optimum compression level. In addition, the anthesis date was incorporated as a cofactor in the analysis. As reported in other studies [29,30,31,32], an adjusted –log10 P > 3 was established as a threshold for considering a marker-trait association (MTA) statistically significant. A moderate threshold at –log10 P > 2.5 was also established for GY. Confidence intervals (CI) for MTAs were calculated for each chromosome according to the linkage disequilibrium (LD) decay reported by Rufo et al. [23]. In order to simplify the MTA information, the associations were grouped into QTL hotspots when at least two MTAs belonging to different traits overlapped their CIs. Circular Manhattan plots were performed using the R package “CMplot” (http://www.r-project.org).

2.6. Gene Annotation

Gene annotation within the CIs of the QTL hotspots was performed using the gene models for high-confidence genes reported for the wheat genome sequence [33], available at https://wheat-urgi.versailles.inra.fr/Seq-Repository/Assemblies. Markers flanking the CIs were used to estimate physical distances from genetic distances.

3. Results

3.1. Phenotypic Data of Root Traits

A summary of the genetic variation of the root traits is shown in

Table 1. Statistics of the seminal root traits.

Table	TRN (N)	RDW (mg)	SRA (°)	SW (g)	SL (cm)
Min	3.2	43	53.1	0.27	14.9
Max	5.4	20.5	125.9	0.63	30.2
Mean	4.4	11.7	98.6	0.48	22.0
SD	0.5	2.2	13.3	0.05	2.8
CV (%)	10.9	18.8	13.5	10.4	12.7
H (%)	52.0	48.5	70.0	75.4	50.0

SD, standard deviation; CV, coefficient of variation; H, repeatability; TRN, total root number; RDW, root dry weight; SRA, seminal root angle; SW, seed weight; SL, seed length.

. The genotypes showed a low coefficient of variation (CV) with a narrow range of variation among traits, from 10.4 for SW to 18.8 for RDW and repeatability (H) ranging from 48.5% for RDW to 75.4% for SW.

The ANOVA (

Table 2. Percentage of the sum of squares of the ANOVA in a set of 55 bread wheat landraces structured into three genetic subpopulations with membership coefficient q > 0.8.

Source of Variation	df	TRN	RDW	SRA	SW	SL
Replicate	1	0.1	0.7	0.0	0.1	0.3
Genotype	54	73.1 ***	64.8 *	82.1 ***	88.8 ***	63.5 *
Between SPs	2	25.3 ***	0.6	10.3 ***	18.9 ***	6.0 *
Within SPs	52	47.8 *	64.2 *	71.8 ***	69.9 ***	57.5 *
WM	16	28.3 ***	31.6	40.8 ***	38.3 ***	26.8 *
NM	14	16.1	14.8	34.4 ***	6.3	25.2
EM	22	55.6	53.6 *	24.8	55.4 ***	48.0
Replicate x Genotype	54	26.8	34.5	17.9	11.1	36.2
Total	108

WM, western Mediterranean; NM, northern Mediterranean; EM, eastern Mediterranean; TRN, total root number; RDW, root dry weight; SRA, seminal root angle; SW, seed weight; SL, seed length. * P < 0.05, *** P < 0.001.

) for cultivars with a high membership coefficient (q > 0.8) showed that for all traits the total variability was mainly explained by the genotype effect, in a range from 63.5% for SL to 88.8% for SW. When the sum of squares of the genotype effect was partitioned into differences between and within SPs, the results revealed that the genetic variability was mainly explained by differences within SPs in a range from 47.8% for TRN to 71.8% for SRA (Table 2). Differences between SPs were statistically significant for SRA, TRN, SW, and SL, in a range from 6.0% of the genotype effect for SL to 25.3% for TRN (Table 2). The sum of squares within SPs was partitioned into western (WM), northern (NM), and eastern (EM) effects, being statistically significant for SRA (40.8%), TRN (28.3%), SW (38.3%), and SL (26.8%) in the western SP, SRA (34.4%) in the northern SP and RDW (53.6%) and SW (55.4%) in the eastern SP.

The ANOVA for grain yield revealed that the genotype effect was the most important in the phenotypic expression of traits, accounting for 59% of the total phenotypic variation, whereas the year effect accounted only for 5%. The interaction accounted for almost 36% of the phenotypic variation although it was not significant (

Table 3. Percentage of the sum of squares for grain yield of the ANOVA in the collection of 170 bread wheat landraces.

Source of Variation	df	Grain Yield	P
Genotype	169	59.2	< 0.001
Year	2	5.1	< 0.001
Genotype x Year	338	35.7	No significant
Total	509

).

The landraces from northern Mediterranean countries showed the highest number of seminal roots with a root angle not statistically different from the western Mediterranean ones. On the other hand, eastern Mediterranean landraces showed the lowest number of roots but the widest angle. These landraces reported the lowest SW and the longest shoots. No differences were reported for RDW among the three SPs (

Table 4. Means comparison of seminal root traits measured in a set of 55 Mediterranean wheat landraces structured into three genetic subpopulations [23] with q > 0.8. Means within columns with different letters are significantly different at P < 0.05 following a Tukey test.

	TRN (N)	RDW (mg)	SRA (°)	SW (g)	SL (cm)
Northern Mediterranean	4.7 a	0.011 a	98.5 b	0.50 a	20.8 b
Western Mediterranean	4.3 b	0.011 a	96.2 b	0.49 a	21.4 ab
Eastern Mediterranean	4.0 c	0.011 a	106.5 a	0.45 b	22.5 a

TRN, total root number; RDW, root dry weight; SRA, seminal root angle; SW, seed weight; SL, seed length.

).

Correlation coefficients between root traits were calculated, showing highly significant correlation coefficients between RDW and SW and RDW and SL (r = 0.47 and 0.45 respectively; P < 0.0001). Moderate significant correlations were reported for TRN with RDW, SW and SRA (r = 0.20, 0.28 and 0.28, respectively), and for SW with SL (r = 0.27). Finally, a negative correlation coefficient (r = −0.12) was found between SRA and SW (Figure 2). GY showed a moderate significant correlation with TRN and SW (r = 0.28 and 0.29, respectively; P < 0.0005).

3.2. Marker-Trait Associations

After filtering for duplicated patterns, missing values, and minor frequency alleles, a total of 10,458 SNPs were used to genotype the panel of 170 wheat landraces [23].

The results of the GWAS for root related traits are reported in Figure 3 and Supplementary Materials, Table S2. Using a common threshold of −log10 P > 3, as reported by other authors [29,30,31,32], a total of 135 MTAs were identified for the analyzed traits. Of these, 50 MTAs corresponded to SW, 39 to RDW, 18 to SL, 17 to SRA, and 11 to TRN. The A and B genomes harbored 46% and 48% of MTAs, respectively, whereas the D genome harbored only 6% of MTAs. The number of MTAs per chromosome ranged from 1 in chromosomes 4D, 5D, and 6D to 14 in chromosome 1B, with a mean of 7 MTAs per chromosome. Most of the MTAs (88%) showed a phenotypic variance explained (PVE) by each MTA in a range of 5% to 10%, and only 2% showed a PVE higher than 15%. Among traits, the PVE mean was stable in a range of 7% (SL) to 9% (RDW).

In order to identify and summarize the genomic regions most involved in trait variation, QTL hotspots were defined when two or more MTAs from different traits were grouped together within the same LD block. LD was previously estimated for locus pairs in each chromosome, and its decay was set to 1 to 10 cM depending on the chromosome [23]. Using this approach, 15 QTL hotspots grouping 43 MTAs were identified (

Table 5. Root QTL hotspots. Positions are indicated in cm.

QTL Hotspot	MTAs	Trait	Chromosome	Peak	CI Left	CI Right
rootQTL1B.1	2	RDW, SW	1B	70.6	69.6	71.6
rootQTL1B.2	2	RDW, TRN	1B	77.5	75.9	79.1
rootQTL1B.3	2	RDW, SL	1B	83.0	81.4	84.6
rootQTL2A.1	2	RDW, SW	2A	47.8	46.7	48.9
rootQTL2A.2	2	RDW, SW	2A	104.1	103.6	104.6
rootQTL2A.3	2	SW, SL	2A	177.5	176.9	178.2
rootQTL2B.1	2	SW, SL	2B	109.5	109.0	110.0
rootQTL3B.1	3	RDW, SW, TRN	3B	62.3	61.4	63.2
rootQTL3B.2	2	SRA, SW	3B	80.6	79.6	81.5
rootQTL5A.1	2	RDW, SL	5A	56.5	56.0	57.0
rootQTL5B.1	2	RDW, SL	5B	95.7	94.5	96.9
rootQTL6A.1	2	RDW, SL	6A	45.8	40.0	51.6
rootQTL6A.2	2	RDW, TRN	6A	76.7	70.7	82.6
rootQTL6A.3	2	RDW, SW	6A	138.4	132.3	144.6
rootQTL7A.1	3	RDW, SRA, TRN	7A	216.6	215.3	218.0

TRN, total root number; RDW, root dry weight; SRA, seminal root angle; SW, seed weight; SL, seed length.

), while 92 MTAs remained as singletons.

The results of the GWAS for GY are reported in Figure 3 and Supplementary Materials, Table S3. A common threshold of −log10 P > 3, detected a total of 40 MTAs, thus a moderate threshold at –log10 P > 2.5 was applied, increasing the number of significant associations to 112. Of these, 32 MTAs corresponded to the year 2016, 30 to 2017, 18 to 2018, and 32 across years. The A and B genomes harbored 43% and 38% of MTAs, respectively, whereas the D genome harbored only 18% of MTAs. The number of MTAs per chromosome ranged from 1 in chromosomes 3D and 6B to 16 in chromosome 1D. Chromosomes 1A, 4D, 5D, and 7D did not show any association. All of MTAs showed a phenotypic variance explained (PVE) by each MTA in a range from 5% to 11%. Most of the MTAs with a PVE > 8% were located on chromosome 1D (76%, 13 out of 17), whereas the percentage increased to 80% among MTAs with a PVE > 10% (4 out of 5).

In order to identify and summarize the genomic regions with a pleiotropic effect for root traits and grain yield, QTL hotspots were defined as previously but including the MTAs for GY. Using this approach, 17 QTL hotspots grouping 81 MTAs were identified (

Table 6. QTL hotspots including grain yield. Positions are indicated in cm.

QTL Hotspot	MTAs	Trait	Chromosome	Peak	CI Left	CI Right
QTL yield/root_1B.1	3	GY, SRA	1B	8.4	7.4	9.4
QTL yield/root_1B.2	3	GY, SW	1B	43.9	42.9	44.9
QTL yield/root_1B.3	3	GY, SW	1B	63.5	61.5	65.5
QTL yield/root_1B.4	8	GY, RDW, SL	1B	83.3	82.3	84.3
QTL yield/root_2A.1	8	GY, RDW, SW	2A	104.1	103.6	104.6
QTL yield/root_2A2	3	GY, SRA	2A	151.3	150.8	151.9
QTL yield/root_3A.1	8	GY, RDW, SL	3A	84.3	81.9	86.7
QTL yield/root_3B.1	5	GY, SW	3B	72.8	70.8	74.8
QTL yield/root_3B.2	3	GY, SRA, SW	3B	80.5	79.6	81.5
QTL yield/root_4B.1	4	GY, SW	4B	76.6	74.1	79.2
QTL yield/root_5B.1	3	GY, SL	5B	57.8	56.8	58.9
QTL yield/root_5B.2	7	GY, RDW	5B	77.3	75.9	78.8
QTL yield/root_5B.3	2	GY, RDW	5B	176.2	175.2	177.2
QTL yield/root_6A.1	5	GY, RDW, SL	6A	45.6	39.6	51.6
QTL yield/root_6A.2	10	GY, RDW, SW	6A	76.6	70.7	82.6
QTL yield/root_7A.1	3	GY, RDW, SW	7A	135	133.2	136.8
QTL yield/root_7B.1	3	GY, RDW	7B	70.0	67.8	72.3

TRN, total root number; RDW, root dry weight; SRA, seminal root angle; SW, seed weight; SL, seed length; GY, grain yield.

). From them, five were in common with those reported only with root traits (rootQTL1B.3, rootQTL2A.2, rootQTL3B.2, rootQTL6A.1, and rootQTL6A.2). GY shared 8 genomic regions with SW and 9 with RDW, 4 with SL, and 3 with SRA, whereas no regions were in common with TRN. In 59% of these genomic regions, GY co-localize with only one root trait, whereas the other 41% co-localize with two different root traits.

In order to identify the most useful markers for selecting for the root traits, extreme phenotypes were identified in the upper and lower 10th percentile of genotypes within the collection for each trait (Figure 4). Among the most significant MTAs for each trait, markers with different alleles between extreme genotypes were identified (

Table 7. Selected significant markers from the GWAS with different allele composition for the upper (UP) and lower (LOW) 10th percentile of genotypes. Different letters on the UP and LOW 10th phenotype indicate that means are significantly different at P < 0.01 following a Tukey test.

Trait	Phenotype			Marker	Chr	Position	R2	Most Frequent Allele
	Mean	UP 10th	LOW 10th					UP	Freq	LOW	Freq
TRN (N)	4.4	5.2a	3.4b	Excalibur_c8613_1266	1B	78.1	0.09	T	0.70	C	0.80
				Tdurum_contig12995_792	5B	98.3	0.12	T	1.00	Y	0.50
RDW (g)	0.012	0.016a	0.007b	RAC875_c8482_160	1A	95.2	0.14	A	0.60	A	0.90
				GENE-0392_97	1A	95.5	0.11	T	0.60	C	0.70
				RAC875_c39401_502	1A	104.1	0.15	T	1.00	C	1.00
				BS00023071_51	1B	95.5	0.15	G	1.00	A	1.00
				Kukri_c14877_303	6A	79.1	0.14	C	0.70	T	1.00
SRA (°)	98.65	119.98a	69.88b	Kukri_c34453_332	4B	70.9	0.11	A	0.60	C	0.70
				BS00028183_51	5B	51.0	0.11	T	1.00	C	0.60
SW (g)	0.487	0.59a	0.37b	RAC875_c37540_583	2B	139.5	0.14	T	1.00	C	0.77
				Ku_c101046_1063	4B	77.7	0.11	Y	0.50	C	0.80
				BS00094770_51	4D	80.4	0.12	C	1.00	A	0.60
				GENE-2220_165	7D	161.1	0.12	T	1.00	T	0.70
SL (cm)	22.04	27.72a	16.80b	Tdurum_contig8158_269	1B	82.4	0.09	A	0.80	G	0.90
				RFL_Contig5277_888	2A	177.6	0.09	C	1.00	C	0.60
				BobWhite_c6966_236	5A	130.9	0.09	T	0.80	C	0.90
				BS00023161_51	5B	48.3	0.08	T	0.90	Y	0.50
				RAC875_rep_c106439_1159	6A	48.1	0.09	A	0.88	C	0.70

Chr, chromosome. TRN, total root number. RDW, root dry weight. SRA, seminal root angle. SW, seed weight. SL, shoot length.

, Figure 5). The frequency of the most common allele among genotypes from the upper 10th percentile ranged from 78% for RDW to 88% for SW, while for the lower 10th percentile it ranged from 65% for TRN and SRA to 92% for RDW (Figure 2).

3.3. Gene Annotation

As reported in Supplementary Materials, Table S4, a total of 1489 gene models were identified within the 15 QTL hotspots using the high-confidence gene annotation from the wheat genome sequence [33]. Genetic distances were converted into physical distances using the position of common flanking markers on the genetic map [24] and the genome sequence. The number of gene models ranged from 224 in rootQTL_2A.2 to 9 in rootQTL_5B.1. Based on the high number of gene models, a selection was made according to gene families involved in root traits, growth and development, and abiotic stress resistance (

Table 8. Selected gene model families.

Description	N QTL Hotspots	Function
F-box family protein	12	Salt and drought stress responses
RING/FYVE/PHD zinc finger protein	12	Salt and drought stress responses
MYB-related transcription factor	8	Salt and drought stress responses
NAC domain-containing proteins	8	Induced by biotic and abiotic stresses
Cytochrome P450 family protein	5	Involved in seed size
BZIP transcription factor	5	Regulated by abiotic stress
Ethylene-responsive transcription factor	4	Induced by biotic and abiotic stresses
Calmodulin	4	Heat shock transduction pathway
Peroxidase	4	Root growth
ABC transporter	4	Control root development
Nucleoside triphosphate hydrolase	3	Associated with drought stress
E3 ubiquitin-protein ligase	3	Associated with drought stress
Glycine-rich protein	2	Enhance drought stress tolerance
Xyloglucan endotransglucosylase/hydrolase	2	Response dehydration, salinity, cold
Aquaporin	2	Drought stress tolerance
Expansin protein	2	Drought tolerance in wheat
Trihelix transcription factor	2	Stomatal development, drought
VQ motif family protein	2	Involved in seed size
Heat shock family protein	2	Induced by abiotic stress
Protein root UVB sensitive 6	2	Early seedling morphogenesis
SAUR-like auxin-responsive family protein	2	Maintain growth during abiotic stress
Bax inhibitor-1 family protein	1	Tolerance to abiotic stresses
Formin-like protein	1	Structure organization in drought-stressed plants
Late embryogenesis abundant protein	1	Participate in drought response
Cell wall invertase	1	Downregulated by drought
Senescence regulator	1	Related to drought stress
Plastid-lipid associated protein PAP/fibrillin	1	Induced by drought
Protein STAY-GREEN LIKE, chloroplastic	1	Improves drought resistance
PI-PLC X domain-containing protein	1	Induced by abiotic stresses
Histidine-containing phosphotransfer protein	1	Enhance tolerance to drought stress
Phospholipase D	1	Enhance drought stress tolerance

). Thus, 31 gene families with a total of 96 gene models remained for subsequent analysis. Among them, F-box and zinc finger family proteins were identified in 12 of the 15 QTL hotspots, whereas 10 gene families were present in only one QTL hotspot. Among chromosomes with QTL hotspots, chromosome 2A had the highest number of gene models (22), whereas chromosomes 5A and 5B had the lowest number (4).

4. Discussion

Breeding for drought adaptation is one of the main challenges to be addressed in the coming years in order to increase wheat production and ensure sufficient food supply in the current scenario of climate change. Roots are crucial in this adaptation, as they are responsible for water and nutrient uptake. The wide morphological plasticity of the root system to different soil conditions and the role of root traits in drought environments are well known [34,35]. Wheat roots reduce their growth in water-limited conditions but increase the water uptake rate, extracting the water from deep soil layers [36]. The shape and spatial arrangement of the RSA can provide a growth advantage and increasing yield performance during periods of water scarcity [37]. Thus, it is necessary to increase the knowledge of the genetics of root architecture in order to improve wheat yield stability under stress conditions by introgressing favorable alleles through breeding programs.

The current study evaluated root-related traits in a collection of Mediterranean bread wheat landraces representative of the variability existing for the species in the Mediterranean Basin [23] with the aim of providing QTL information for these traits regarding seminal roots. Seminal roots are important for early vigor and crop establishment in dryland areas because they explore the soil for nutrients and water [38]. Moreover, it has been reported that under drought stress, seminal roots activity is more important than that of nodal roots [39]. Additionally, field phenotyping of hundreds of genotypes is a complex and expensive task. As the root geometry of adult plants is strongly related to the SRA [5], it may be assumed that genotypes that differ in root architecture at an early developmental stage would also differ in the field at later growth stages, when nutrient and/or water capture become critical for yield performance [8].

The range of variation for the traits analyzed in the present study (from 10.9% for TRN to 18.8% for RDW) is in agreement with those reported for elite durum wheat cultivars by Canè et al. [8], who explained this variability as an adaptive value for the environmental conditions of the region of origin of the cultivars. Moreover, the high repeatability found for the traits supports the approach followed to analyze the seminal roots under controlled conditions.

Landraces from the eastern Mediterranean Basin showed the widest SRA, the lowest SW, the longest SL, and the lowest number of roots. According to previous studies in durum wheat [18,40], landraces from southeastern Mediterranean countries corresponding to the warmest and driest areas of the Mediterranean Basin, reported more grains per unit area and lighter grains than those developed in cooler and wetter zones of the region. Although it has been reported that in water-limited environments a vigorous root system could have benefits at the beginning of the growing season because it offers a more efficient water capture [41], no significant differences were observed for RDW among the SPs in the current study. Moreover, our results for SRA are in agreement with those reported by Roselló et al. [18], who found that durum wheat landraces from the eastern Mediterranean have the widest root angle, which probably allows them to cover a larger soil area and be more efficient in water uptake than landraces that originated in wetter areas.

Although not significant, probably due to the very early stage when the root traits were measured, the correlation between SRA and SW was negative. The same result was also reported by Canè et al. [8], who suggested that it could be due to the influence of the root angle on the distribution of the roots on soil layers and, therefore, the water uptake from deeper layers. On the other hand, the correlation between RDW and SW was positive, in agreement with the findings of Fang et al. [42], thus indicating the effectiveness of greater root mass for obtaining more soil water for plant growth and grain filling in drought. Seedling growth has also been related to SW in wheat [43]. The vertical distribution of the root system can have a strong effect on yield [44], so mass root concentrated in upper layers can be more effective for resource capture, while roots in deeper layers have more access to deep water.

The complexity of the genetic control of root traits was confirmed with 135 marker-trait associations identified in the current study. Their distribution across genomes was similar in the A and B genomes (46% and 48%, respectively), leaving only 6% of MTAs in the D genome. These results agree with the lower genetic diversity and higher LD found in the D genome, as reported previously [23]. According to Chao et al. [45], the different levels of diversity in wheat genomes could be due to different rates of gene flow from the ancestors of wheat, since polyploidy bottleneck resulting from speciation reduced diversity and increased the levels of LD in the D genome in comparison with the A and B genomes.

In order to simplify and to integrate closely linked MTAs in a consensus region, QTL hotspots were identified based on the results of LD decay reported in [23]. LD decay was used to define the CIs for the QTL hotspots. Following this approach, 43 MTAs were grouped in 15 QTL hotspots. The genomic position of QTL hotspots was compared with previous studies reporting meta-QTLs for root traits [46] and MTAs from GWAS studies in order to detect previously identified regions controlling root traits. Among the 15 QTL hotspots, only rootQTL6A.3 was located in the same region of a previously mapped meta-QTL, RootMQTL74 [46]. When compared with MTA-QTLs reported by [18] in durum wheat Mediterranean landraces, the QTL hotspot rootQTL6A.3 corresponded to the MTA-QTLs mtaq-6A.3 and mtaq-6A.6. This hotspot was also in the same region of a major SRA QTL identified by Alahmad et al. [47] and by a QTL controlling root growth angle identified by Maccaferri et al. [9], who also found a QTL for grain weight that is located in a common region with the hotspot rootQTL2A.2, which includes an MTA for SW. rootQTL3B.1 shared a common position with an MTA reported by Ayalew et al. [48] on chromosome 3B under stress conditions. rootQTL7A.1, including an MTA for RDW, was located in a similar position as MLM-RDWB-10 reported by Li et al. [49] and associated with RDW at the booting stage. Finally, no genomic regions were shared with the study carried out by Beyer et al. [50]. Only four of the 15 QTL hotspots identified in this work had been detected previously, suggesting the importance of wheat Mediterranean landraces for the identification of new loci controlling root-related traits.

As reported in previous studies, at early developmental stages [8,18] the co-location of MTAs for grain yield and root related traits within the same QTL hotspot suggests their pleiotropic effect, however, deeper analyses should be necessary to confirm it. In durum wheat elite cultivars, Canè et al. [8] found that 30% of the QTLs affecting root system architecture were included within QTLs for agronomic traits. More recently, Roselló et al. [18] using a collection of Mediterranean durum wheat landraces found that 45% of QTL hotspots for root related traits were mapped in similar regions to yield-related traits reported for the same collection of landraces.

From a breeding standpoint, exploiting genetic diversity from local landraces is a valuable approach for recovering and broadening allelic variation for traits of interest [19]. Therefore, identifying the genotypes showing the extreme phenotypes within the pool of Mediterranean landraces and the associated markers provide the opportunity for introgressing suitable traits in elite cultivars by marker-assisted breeding using the most recent technologies to speed the process.

The availability of a high-quality reference wheat genome sequence [33] enabled us to quickly identify gene models corresponding to QTLs. Thus, the genetic position of the CIs of the QTL hotspots was projected into physical distances on the reference sequence to search for putative candidate gene models. To narrow the number of candidates, only gene models involved in the development and abiotic stress according to the literature were taken into consideration. Therefore, of 1489 gene models identified within the 15 QTL hotspots, only 31 gene families were selected.

F-box and zinc finger family proteins were the most represented, each one appearing in 12 hotspots. F-box proteins play important roles in plant development and abiotic stress responses via the ubiquitin pathway [51] and the ABA signaling pathway [52]. In wheat, the F-box protein TaFBA1 is involved in plant hormone signaling and response to abiotic stresses and is expressed in all plant organs, including roots [53]. The overexpression of TaFBA1 in transgenic tobacco reported by Li et al. [54] to improve heat tolerance resulted in increased root length in the transgenic plants. Zinc finger proteins are involved in several processes, such as regulation of plant growth and development, and response to abiotic stresses [46]. In Arabidopsis and rice, they play a role in tolerance to drought and salt stresses [55], while in wheat the overexpression of TaZFP34 enhances root-to-shoot ratio during plant adaptation to drying soil [56].

Other kinds of gene models found in a high number of QTL hotspots were MYB transcription factors and NAC domain-containing proteins, each of them presents in 8 hotspots. MYB domain-containing transcription factors are involved in salt and drought stress adaptation in wheat. Some examples in wheat are the genes TaMyb1, TaMYBsdu1, and TaMYB33. The expression of TaMyb1 in roots is strongly related to responses to abiotic stresses [57]. The gene TaMYBsdu1 was found to be upregulated in leaves and roots of wheat under long-term drought stress [58]. Finally, the overexpression of TaMYB33 in Arabidopsis enhances tolerance to drought and salt stresses [59]. NAC domain-containing proteins have been described to play many important roles in abiotic stress adaptation [46]. Xie et al. [60] reported that NAC1 promoted the development of lateral roots. Similarly, He et al. [61] found that the expression of AtNAC2 in response to salt stress led to an increase in the development of lateral roots. Xia et al. [62] demonstrated that the gene TaNAC4 is a transcriptional activator involved in wheat’s response to biotic and abiotic stresses.

Proteins belonging to the cytochrome P450 family and bZIP transcription factors were present in five QTL hotspots. The first class of proteins belongs to one of the largest families of plant proteins, with genes affecting important traits for crop improvement such as TaCYP78A3, which is involved in the control of seed size [63]. bZIP transcription factors are involved in abiotic stress response [64]. In Arabidopsis, it has been observed that the overexpression of TabZIP14-B, involved in salt and freezing tolerance, hindered root growth in transgenic plants in comparison with the control plants [65].

Other proteins involved in root growth and development are the peroxidases and ABC transporters that were identified in four QTL hotspots. Extracellular peroxidases are involved in plant defense reactions against biotic and abiotic stresses through the generation of reactive oxygen species in wounded root cells [66]. In Arabidopsis, the ABC transporter AtPGP4 is expressed mainly during early root development, and its loss of function enhances lateral root initiation and root hair development [67]. Gaedeke et al. [68] reported a new member of the ABC transporter superfamily of Arabidopsis thaliana, AtMRP5. Using reverse genetics, these authors found that the recessive allele mrp5 exhibited decreased root growth and increased lateral root formation. In addition to peroxidases and ABC transporters, other proteins identified in four QTLs were the ethylene-responsive transcription factors (ERFs), found to be involved in the response to abiotic stresses. In wheat, the ERF TaERFL1a is induced in wheat seedlings in response to salt, cold, and water deficiency [69].

Other family proteins involved in drought stress, seed size, or early development were represented in a lower number of QTL hotspots. Among them, aquaporins are known to affect drought tolerance influencing the capacity of roots to take up the soil water [70]. The expansins were suggested to be involved in root development, as the overexpression of the wheat expansin TaEXPB23 improved drought tolerance by stimulating the growth of the root system in tobacco [71].

5. Conclusions

The exploitation of unexplored genetic variation present in local landraces can potentially contribute to breeding programs aimed at enhancing drought tolerance in wheat. Roots are crucial for adaptation to drought stress because they are the plant organ responsible for water and nutrient uptake and interaction with soil microbes. Thus, designing and developing novel root system ideotypes could be one of the targets of wheat breeding for the coming years. The variability found in the Mediterranean wheat landraces together with the newly identified QTL hotspots shows landraces as a valuable source of favorable root traits to introgress into adapted phenotypes through marker-assisted breeding. Among the different marker trait associations, those reported in extreme genotypes could result as a starting point to develop new mapping populations to fine map the corresponding traits.