Ten agronomic and quality traits were assessed for three durum wheat elite line panels: days to heading (days, DTH); plant height (cm, PH); lodging (%, LOD); yellow rust (%, YR); yellow colour (YC); sedimentation index (cm3, SDS); sedimentation volume (cm3, SV); grain protein content (%, GPC); thousand kernel weight (g, TKW); and grain yield (Kg/ha, GY). Agronomic traits (DTH, PH, LOD, YR, TKW and GY) were assessed under both water treatments (FI and RI), and quality traits (YC, SV, SDS and GPC) were only evaluated under FI conditions. Visual disease evaluation and phenology assessments were made in the field, while quality parameters were evaluated on grain samples post-harvest. DTH, PH and LOD and YR were visually assessed at the field trials in Yaqui, while YR was assessed at Toluca. To assess DTH, heading date was recorded as the time when 50% of the spikes have emerged from the flag leaf sheath (stage 59 in Zadoks scale [51]); PH was recorded by measuring the distance between the stem's base and the top of the spike (awns not included); LOD was assessed as the percentage of lodging plot; and YR was assessed as the percentage of leaves with rust pustules. YC was assessed by a rapid colorimetric method with a Minolta color meter following CEN/TS 15,465:2008 [52–54]; SDS was evaluated following UNE 34,903:2014 [55,56]; SV and GPC were assessed by Near-infrared spectroscopy (NIRs) [57]; TKW was measured by weighing 2 samples of 100 entire kernels randomly chosen previously dried at 70 ◦C for 48 h.

The correlation between the assessed traits was analysed using the 'cor' function in R [58–60]. Then, an analysis of variance (ANOVA) was undertaken, using the 'aov' function in R [61], to obtain the descriptive statistics for each trait.

Traits were analysed using a Q + K linear mixed-model [62,63] which follows the model equation:

$$\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \mathbf{S}\boldsymbol{\alpha} + \mathbf{Q}\mathbf{v} + \mathbf{Z}\boldsymbol{\mu} + \boldsymbol{\varepsilon} \tag{1}$$

where y is a vector of observed phenotypes; X, S and Z are matrices related to β, α and μ, respectively; β is a vector of fixed effects; α is a vector of marker effects; Qv is a vector of population effect; μ is a vector of polygenic effects (with covariance proportional to a kindship or relationship matrix); and ε is a vector of residuals.

These analyses were carried out using GenStat (14th Edition) to generate the best linear unbiased estimates (BLUEs) of variety performance in different ways: (i) across years and blocks; (ii) across years for each block (FI and RI); (iii) across a reduced dataset (years 2013 and 2014) and blocks; and (iv) across the reduced dataset for each block. The resulting datasets (available in Supplementary Material Table S2) were then used in different association mapping analyses.

#### *2.3. Population Structure and Linkage Disequilibrium*

Population structure was assessed using principal component analysis (PCA) based on the combined DArT and SNP genotyping datasets. Euclidean distances were calculated using the R package 'ggfortify' [64] and the PCA was visualised with 'ggplot2' [65].

The pattern of linkage disequilibrium (LD) was assessed between each pair of SNP markers on the same chromosome across the two constitutive genomes with the allele frequency correlation (r2) using the 'popgen' package in R [66]. A heatmap was obtained with the *D'* and r2 values for each chromosome and a scatterplot to determine LD decay (genetic distance in cM).

#### *2.4. Association Mapping (AM)*

The AM analyses were performed on the BLUEs obtained above using an additive model with 'rrBLUP' [67] and 'GWASpoly' [68] packages in R in different ways. Two marker-based kinship matrices (k-matrix), created from a subset of 14,588 DArT and 5716 SNP markers, respectively, were used for the adjustment based on relatedness of individuals (Supplementary Materials Tables S3 and S4). A minor allele frequency (maf) threshold of 0.05 was used. To establish a *p*-value detection threshold

for statistical significance of associations, the Bonferroni correction, which employs a threshold of α/m to ensure the genome-wide type error I of 0.05, was applied with a total of 1000 permutations.

Associated DartSeqTM and SNP markers were blasted against the wheat reference assembly RefSeqv1 [69] with no indels or mismatches allowed, using an ad hoc Java program, to confirm their physical mapping location on the A or B genomes. The molecular markers were also mapped against the durum wheat genome (https://www.interomics.eu/durum-wheat-genome) to confirm their physical positions. In addition, to identify candidate genes, results were filtered, selecting the hits with best *e*-value for each molecular marker, and candidate genes were manually selected based on their annotations.

#### **3. Results**

#### *3.1. Phenotypic Assessment*

Results from the ANOVA for all the traits across years and water treatments are shown in Table 3. The mean phenotypic values across years were calculated for each block and panel to evaluate the influence of water conditions on the assessed traits (Table 3). Days to heading during the field seasons assessed ranged from 63 to 94 days. In plots with lower water availability (RI block), the spike emergence from the flag leaf took place approximately 11 days earlier than in FI plots. Plant height ranged from 39 to 110 cm showing differences between water regimes, with a decrease of 25–30 cm under RI conditions. Likewise, and as result of the RI treatment, GY (ranging from 1.35 to 10.63 ton/ha) and TKW (from 29.6 to 63.2 g) also varied, being reduced by 4–5 tons/ha and 7–10 g, respectively, in the low water availability RI treatment. This strong RI treatment resulted in very low heritability values for DTH, PH and LOD.

Several significant phenotypic correlations were observed between the analysed traits (Figure 1 and Supplementary Materials Table S5). The most correlated traits were PH and GY (r = 0.90, *p*-value = <2.2 × 10−16), followed by DTH and GY (r = 0.87, *p*-value = <2.2 × 10−16), SDS and SV (r = 0.85, *p*-value = <2.2 × 10−16) and also DTH and PH (r = 0.82, *p*-value = <2.2 × 10−16). Other traits showed important positive correlations too, including YC and DTH (r = 0.69, *p*-value = 2.82 × 10−09), GY and TKW (r = 0.66, *p*-value = <2.2 × 10−16), PH and TKW (r = 0.62, *p*-value = <2.2 × 10−16), GY and YC (r = 0.53, *p*-value = 0.039) and DTH and TKW (r = 0.49, *p*-value = <2.2 × 10−16). Negative correlations were also observed for SDS and GPC (r = −0.38, *p*-value = <2.2 × 10−16), and for GPC and YC (r = −0.08, *p*-value = <2.2 × 10−16) (Figure 1).

**Figure 1.** Phenotypic correlations. YR: yellow rust; DTH: days to heading; PH: plant height; GY: grain yield; TKW: thousand kernel weight; YC: yellow colour; SV: sedimentation volume; SDS: sedimentation index; and GPC: grain protein content.

#### *3.2. Population Structure and Linkage Disequilibrium*

The PCA used a total of 14,588 DArT and 5716 SNP markers. The first and second principal components explained 3.91% of the genetic variation (Figure 2). No underlying genetic structure was detected within or between the panels assessed. LD was estimated using the mapped SNP markers dataset. LD decay was determined within 20–30 cM for all the chromosomes (Figure 3). Using the classification defined by Maccaferri et al. [70], the markers presented loose linkage (Class 2), showing a distance value between 21 to 50 cM.

**Figure 2.** Principal components analysis (PCA) of the genotypic data.

The LD decay was established

 between 20–30 cM.

#### *3.3. AM Analysis*

Thirty-seven significant marker-trait associations (MTAs) were detected for TKW and SV across all years and water treatments with most of the significant markers located on chromosome 2A (Table 4). Twenty DArT and seven SNP markers were found in association with TKW on chromosome 2A (with additive effects ranging from −3.41 to 3.46). In addition, eight unmapped DArT and one SNP marker were also associated with TKW (additive effects ranged from −3.39 to 3.46 g). Most of these MTAs showed a negative additive effect, reducing the final weight value (ranging from −2.84 to −3.19 g), and only two MTAs were found to increase TKW (values of 2.97 and 3.09 g). Finally, a single SNP associated with SV was located on chromosome 1B (showing a positive effect increasing the final value by 1.26 mL). The resulting Manhattan and QQ-plots from this AM analysis are included in Supplementary Materials Figures S2–S5.

The AM analyses on partitioned subsets of the data consistently detected the QTLs for TKW and SV. Nevertheless, the individual assessment of the water treatments significantly reduced the number of MTAs found, due in part to less available data for the RI block (Supplementary Material Table S6). The initial dataset of 294 durum wheat elite lines was reduced to 200 lines (assessed during the 2013 and 2014 seasons) to give a dataset balanced across assessment years. Using this reduced dataset for AM analysis, the results confirmed the QTLs previously found for the full dataset (Supplementary Materials Table S6). The analysis also detected an additional locus for SV on chromosome 4A, and a locus for YR on chromosome 1B (with additive effects of −0.84 and 2.79, respectively).

#### *3.4. Candidate Genes Analysis*

The marker *SNP620*, located on chromosome 1B and detected in association with SV, was found included into a cluster of 12 genes encoding blue copper proteins (BCP), with homoeologs in the two durum wheat subgenomes (Figure 4, Table 5 and Supplementary Material Figure S6). In the hexaploid wheat genome, this set of genes form a cluster of homeolog triads [72] with a total of 31 genes (Supplementary Materials Table S7 and Figure S7). Additionally, another interesting gene (*TraesCS1B01G568400LC.1*) was found closer this marker, coding for the ubiquinone biosynthesis *O*-methyltransferase.

There were several markers located in chromosome 2A, in close proximity to some interesting genes. Markers *SNP1183*, *SNP1184* and *DArT3165* were found next to several genes encoding reductase-1 (Figure 4 and Table 5). In addition, the marker *SNP8395*, also located on the same chromosome, was found in proximity to the gene *TraesCS2A01G309700.1*, which encodes a type A response regulator 1 (Figure 4 and Table 5).

Significant MTAs from the partitioned analysis also allowed the identification of potentially interesting genes. On chromosome 1B, marker *DArT1744*, previously described by Mérida-García et al. [73] related to high molecular-weight glutenin subunits, was found in proximity to genes encoding isocitrate dehydrogenase kinase/phosphatase G and leucine-rich repeat receptor-like protein kinase family protein. These genes participate in the carbohydrate metabolism during the Krebs cycle and play a crucial role in plant development and stress responses, respectively [74,75]. On this chromosome, another marker (*SNP809*) was found in proximity to some interesting genes encoding sugar transporter proteins. Additionally, some markers located on chromosome 2A were found in proximity to Acyl-CoA *N*-acyltransferase genes (*SNP1206*, *SNP8395* and *DArT3180*) and chloroplastic zeaxanthin epoxidase (*SNP1189*) (Supplementary Materials Table S8).

**Figure 4.** Location of candidate genes on chromosomes 1B and 2A. Markers are indicated with the symbol '-'; blue copper proteins are shown in blue colour; ubiquinone biosynthesis *O*-methyltransferase in purple colour; the regulator response gene in brown colour; and for reductase 1 genes in green colour.

**Table 5.** Candidate genes for markers with significant marker-traits associations. Genes located in proximity of markers found in association with TKW andacross years and water treatments (within a ±50 kb window). Values for physical position and distance are indicated in base pairs (bp). Chr: chromosome position.Blue copper proteins are shown in blue colour; ubiquinone biosynthesis *O*-methyltransferase in purple colour; the regulator response gene in brown colour; and

 SV  for

#### **4. Discussion**

The maintenance of crop production is a current and pressing need given growing populations and reduced availability of arable land [76]. There is an increasing need for breeding programs to improve yield potential and the adaptation of new varieties to different biotic and abiotic stresses [77]. Abiotic stresses, including drought and heat, are impacting productivity on the million hectares of wheat grown worldwide each year [78]. Detailed molecular and phenotypic characterization are valuable tools in the dissection of complex traits [79], and especially those that are influenced by water availability [14].

The improvement of key traits is essential for better end-use production quantity and quality in wheat [80]. In this study, we analysed a set of 10 agronomic and quality traits under full irrigation conditions (FI), with an additional six traits also assessed under low water availability (RI) in order to understand trait variation under contrasting water regimes in the CIMMYT durum wheat breeding programme. Irrigation conditions influenced some important yield and yield-related traits such as GY and TKW, as well as adaptive traits including DTH and PH (Table 3). The RI treatments had decreased final yields in line with previous observations [6,81]. Previous reports have also shown TKW to be reduced by high temperatures [17], most likely related to water availability. Groos et al. [82] assessed the genetic basis of grain yield and protein content in a segregating population of wheat RILs grown at six locations and also identified effects from GxE interactions involving protein content and yield. Our mean trait values corroborated this, with the highest values for GY recorded for FI blocks across panels (see Table 3). A similar trend was shown for DTH, PH and TKW, which decreased under low-water regimes.

Correlations between the assessed traits showed that GY was positively correlated with two different phenology traits (PH and DTH). This is in agreement with Maccaferri et al. [6], who showed important positive and negative correlations for GY and DTH, and also positive correlations for GY and PH in several environments with different water regimes. DTH and PH were also positively correlated (Supplementary Materials Table S5), with taller plants having a longer time period to the emergence of the tip of the spike stage.

Wheat TKW is an important yield component with a direct effect on grain yield [83,84]. In line with this, our results showed a significant and positive correlation between TKW and GY. However, the previously reported negative correlation between TKW and DTH [6] was not observed, potentially as result of temperatures and water availability from emergence to heading, and also from heading to harvest. Rharrabti et al. [17] previously highlighted a positive correlation between protein content and TKW, which is in agreement with the results obtained in the present study. They highlighted that warm temperatures during grain filling could negatively affect this correlation.

Significant associations between endosperm proteins (gliadin and glutenin subunits) and SV have been previously highlighted [85,86]. Here we found a positive correlation (r = 0.15) between SV and GPC. This correlation is thought to be the result of grain protein content influencing the sedimentation volume value [87], which depends on the degree of protein hydration and oxidation [88]. Finally, for sedimentation index (SDS) analysis, we observed a negative correlation with protein content, in agreement with results presented by Rharrabti et al. [17,45]. This is also in agreement with Oelofse et al. [89] who highlighted the significant influence of protein content on the SDS sedimentation test [90–92].

The SNP and DArT markers used to analyze patterns of genetic structure (Figure 2) and LD (Figure 3) for the durum wheat lines revealed no detectable genetic structure and consistent patterns of LD across chromosomes (LD was estimated to extend ~20 cM). These results support the suitability of durum elite line sets currently in use in breeding programmes for the practical application of GWAs analysis. The rate of unmapped markers was lower for SNP than for DArT markers (27.5% vs. 42.0%, Table 2), indicating higher precision in genetically mapping SNP markers, probably as a result of co-dominance.

In the assessment of MTAs for quality and yield-related traits, different AM analyses were performed on subsets of the dataset. Several MTAs for SV and TKW were detected across years and water regimes, located on chromosomes 1B and 2A, respectively. All GWAS analyses corroborated the major QTLs previously detected, and reported two new QTLs, one for YR in chromosome 1B, and another for SV in chromosome 4A.

Associations on chromosome 1B were significant for wheat quality. There are known loci including Gli-B1/Glu-B3 on this chromosome, which are of great importance for some gliadin and glutenin subunits [93]. In fact, the candidate gene analysis reported the presence of a high molecular-weight glutenin subunit (HMW-GS) gene in the proximity of marker *DArT1744* (found to be significantly associated with SV and SDS), which was previously described by Mérida-García et al. [73] in association with specific weight. In line with this, Pogna et al. [93] highlighted the importance of Glu-B3 for determining protein quality with these endosperm proteins showing significant effects on SV, which also showed a high and positive correlation with SDS in our study (Figure 1 and Supplementary Materials Table S5). Likewise, Blanco et al. [86] reported three QTLs on chromosomes 1B, 6A and 7B (based on the analysis of 259 polymorphic markers) associated with SDS and SV in a recombinant inbred line population. In the present study, we found a SNP marker (*SNP620*) associated with SV, showing a positive additive effect of 1.26 (see Table 4) and also with SDS (marker effect of 0.11) (Supplementary Materials Table S6). This marker was previously placed on chromosome 1B, in the same location as MTAs for gluten index and sedimentation index [73]. Other previous studies, such as Reif et al. [94] and Fiedler et al. [95], also reported markers associated with SV on chromosome 1B, but with differing genetic positions. The additional locus for SV found on chromosome 4A (marker *DArT9459*) has not been previously reported.

Marker *SNP620*, significantly associated with SV, is located within a cluster of homoeolog gene triads coding for blue copper proteins (Table 5 and Supplementary Materials Figure S6). These proteins have been described containing a copper atom, and participate in redox processes [96], with a crucial role in the electron shuttle in plants [97]. In addition, Yao et al. [98] described the blue copper protein genes as the targets of miR408 in wheat, which is involved in the regulation of gene transcription required for heading time [99]. In our study *SNP620* was also found in proximity to a gene coding for an ubiquinone biosynthesis O-methyltransferase. Liu et al. [100] highlighted its crucial role as an electron transporter in the electron transport chain of the aerobic respiratory chain. This ubiquinone gene is involved in plant growth and development, is implied in chemical compounds biosynthesis and metabolism which are involved in plant responses to stress, and also participates in gene expression regulation and cell signal transduction [100].

On chromosome 1B we also found a significant MTA for yellow rust, in agreement with previous studies in durum and bread wheat, which placed different markers significantly associated with this trait, but in differing genetic positions [101–104]. The candidate gene analysis revealed the proximity of this marker (*SNP809*) to sugar transporter protein genes. Sugars are formed during the photosynthetic process and are essential for plant nutrition. The sucrose transport has been considered of great importance for plant productivity [105]. In line with this, the sucrose is involved in the gene expression regulation of the supposedly sugar-sensing pathway [106,107].

The majority of MTAs for TKW were located on chromosome 2A, showing both positive and negative effects. Previous studies have reported different markers in association with this quality trait, including a number mapped on chromosome 2A [38,108–110]. One of the markers found by Yao et al. [38] (SSR marker *gwm445* on chromosome 2A (68.2 cM), belongs to the same QTL found in this study for the marker *SNP1153* (chromosome 2A, 68.6 cM), and also found by Juliana et al. [111] in bread wheat lines from CIMMYT's first year-yield trials. Sukumaran et al. [43] analysed a durum wheat panel of 208 lines under yield potential, heat and drought stress conditions, and identified markers on chromosome 2A with a similar position to those detected in this study (4 markers at 70 cM and 6 markers at 69 cM) under heat stress conditions. They highlighted that several SNP markers were related to transmembranes or were uncharacterized proteins. We found several candidate genes for this

important TKW QTL (Table 5) among which the most striking feature is the presence of four reductase 1 genes (NADPH-dependent 6 -deoxychalcone synthase) and a type A response regulator 1 (Figure 4). These genes are both related with photosynthesis. Hu et al. [112] highlighted that NADPH plays a crucial role in biological processes in plants, such as the regulation of the production of reactive oxygen species (ROS) for the stress tolerance [113,114]. Additional GWAS analyses using reduced datasets revealed other interesting genes for this QTL (chromosome 2A, Supplementary Materials Table S8), encoding for the Acyl-CoA *N*-acyltransferase and the chloroplastic zeaxanthin epoxidase. The first gene has several functions in signaling and metabolic pathways [115]. The zeaxanthin epoxidase plays an important role in the xanthophyll cycle and abscisic acid (ABA) biosynthesis. The xanthophyll cycle has a main function in the dissipation of light energy excess and also increasing the photosynthetic system stability [116].

The proposed approach has successfully detected genetic markers with significant associations with TKW, SV, SDS and YR. These are of potential use in durum wheat breeding programs, and can be further interrogated to the candidate gene level using the RefSeqv1 bread wheat genome reference [69] and the durum wheat genome reference [71].

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4395/10/1/144/s1: Figure S1: Rainfall data for Ciudad Obregon (Mexico) for the growing seasons 2012–2015; Figure S2: Manhattan plots for durum wheat mapped DArT markers. DTH: days to heading; PH: plant height; GY: grain yield; TKW: thousand kernel weight; YC: yellow color; SV: sedimentation volume; SDS: sedimentation index; and GPC: grain protein content; Figure S3: Manhattan plots for durum wheat mapped SNP markers. DTH: days to heading; PH: plant height; GY: grain yield; TKW: thousand kernel weight; YC: yellow color; SV: sedimentation volume; SDS: sedimentation index; and GPC: grain protein content; Figure S4: Quantile quantile-plots from genome-wide association studies (GWAS) analysis for durum wheat DArT markers (mapped and unmapped). DTH: days to heading; PH: plant height; LOD: lodging; GY: grain yield; TKW: thousand kernel weight; YC: yellow color; SV: sedimentation volume; SDS: sedimentation index; and GPC: grain protein content; Figure S5: Quantile quantile-plots from GWAS analysis for durum wheat SNP markers (mapped and unmapped). DTH: days to heading; PH: plant height; LOD: lodging; GY: grain yield; TKW: thousand kernel weight; YC: yellow color; SV: sedimentation volume; SDS: sedimentation index; and GPC: grain protein content; Figure S6: Blue copper protein gene cluster on durum wheat chromosome 1B. High confidence genes are shown in green colour, low confidence genes are shown in yellow; Figure S7: Cluster tree of blue copper protein gene homoeologs in bread wheat (RefSeqv1 [69]). For chromosome 1A, high confidence (HC) and low confidence (LC) genes are shown in brown and orange colour, respectively; for chromosome 1B, HC and LC genes are shown in dark and light green colour, respectively; for chromosome 1D, HC and LC genes are shown in dark and light blue colour, respectively; Table S1: Durum wheat elite lines assessed; Table S2: Best Linear Unbiased Estimates (BLUEs) outputs for all assessed traits and the association mapping analyses performed in durum wheat: [i] across years and blocks; [ii] across years for each block (FI and RI); [iii] across years and blocks for a reduced dataset (years 2013 and 2014); and [iv] across the reduced dataset for each block. DTH: days to heading; GPC: grain protein content; GY: grain yield; PH: plant height; SDS: sedimentation index; SV: sedimentation volume; TKW: thousand kernel weight; and YC: yellow colour; YR: yellow rust; LOD: lodging; Table S3: Kinship matrix for durum wheat DArT markers; Table S4: Kinship matrix for durum wheat SNP markers; Table S5: Phenotypic correlations between the assessed traits in durum wheat and their corresponding p values. YR: yellow rust; DTH: days to heading; PH: plant height; LOD: lodging; GY: grain yield; TKW: thousand kernel weight; YC: yellow color; SV: sedimentation volume; SDS: sedimentation index; and GPC: grain protein content; Table S6: Marker-trait associations found for the association mapping analyses performed in durum wheat: [i] across years and blocks; [ii] across years for each block (FI and RI); [iii] across years and blocks for a reduced dataset (years 2013 and 2014); and [iv] across the reduced dataset for each block. SV: sedimentation volume; TKW: thousand kernel weight; SDS: sedimentation index; YR: yellow rust; "-": unmapped marker; Table S7: Homoeolog triads for blue copper protein genes mapped in the wheat reference assembly RefSeqv1 [69]; Table S8: Candidate genes for GWAS analyses performed in durum wheat: [i] across years for each block (FI and RI); [ii] across years and blocks for a reduced dataset (years 2013 and 2014); and [iii] across the reduced dataset for each block. Molecular markers mapping positions are shown both in the durum wheat genome [71] and the wheat reference assembly RefSeqv1 [69]; Supplementary Material S1. R script used to perform the GWAS analysis; Supplementary Material S2. R script used to perform the LD analysis.

**Author Contributions:** P.H., I.S. and K.A. conceived the experiment. K.A. and I.S. selected the plant materials and agronomic traits. K.A. managed the field trials and contributed the phenotypic data. I.S., P.H., G.D., S.G., R.M.-G. and A.R.B. analysed the genotypic and phenotypic data. R.M.-G. isolated the DNA. R.M.-G., A.R.B. and P.H. carried out the AM analyses. R.M.-G., A.R.B. and P.H. drafted the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by project P12-AGR-0482 from Junta de Andalucía (Andalusian Regional Government), Spain (Co-funded by FEDER). ARB is supported by the UK Biotechnology and Biological Sciences Research Council (BBSRC) 'Designing Future Wheat' cross-institute strategic programme (BB/P016855/1). PH is supported by the Spanish Ministry of Economy, Industry and Competitiveness (MINECO) project AGL2016-77149-C2-1-P.

**Conflicts of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

#### **References**

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*

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

#### **Rubén Rufo 1, Silvio Salvi 2, Conxita Royo 1 and Jose Miguel Soriano 1,\***

Received: 3 April 2020; Accepted: 23 April 2020; Published: 25 April 2020

**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.

**Keywords:** drought stress; association mapping; root system architecture; QTL hotspot; seminal root

#### **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–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–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.

**Figure 1.** Experimental setup for the analysis of seminal root traits. Seeds were placed 8 cm apart on moist filter paper (**A**) and kept in a box with distilled water in a growth chamber for 14 days at 22 ◦C under a 16-h light photoperiod (**B**). (**C**) Example of seminal root angle measurement, using ImageJ software.

#### *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 m2. The experimental design is shown in Supplementary Materials, Figure S1. Sowing density was adjusted to 250 germinable seeds m2. Weeds and diseases were controlled following standard practices at the site. The anthesis date was determined in each plot. Grain yield (GY, t ha−1) 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:

$$\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \mathbf{Z}\boldsymbol{\gamma} + \boldsymbol{\varepsilon} \tag{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–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. 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.

**Table 1.** Statistics of the seminal root traits.

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 ANOVA (Table 2) 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.

**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.

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.

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).

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

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).

**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, 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).

**Figure 2.** Correlations between seminal root traits and grain yield. On the right side are shown the values of the correlation coefficients (*r*). SL, seed length; RDW, root dry weight; SW, seed weight; TRN, total root number; SRA, seminal root angle; GY, grain yield.

#### *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–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).

**Figure 3.** GWAS for root related traits (left circle) and grain yield for 3 years and across years (right circle). From the inside out, root traits correspond to RDW, SW, TRN, SRA, and SL, whereas for GY corresponds to 2016, 2017, 2018 harvesting seasons and the mean across years.

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), 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). 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.

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

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

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

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, 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).

**Figure 4.** Extreme phenotypes for SRA and TRN. The means correspond for 3 observational units of the genotype for SRA and 5 observational units of the genotype for TRN.

**Figure 5.** Marker allele frequency means from landraces within the upper and lower 10th percentile for the analyzed traits. All significant markers shown in Table 5 are included. TRN, total root number; RDW, root dry weight; SRA, seminal root angle; SW, seed weight; SL, seed length.

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

#### *Agronomy* **2020**, *10*, 613

#### *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). 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.

**Supplementary Materials:** Supplementary materials can be found at http://www.mdpi.com/2073-4395/10/5/613/s1. Figure S1: Field scheme of the experimental design for grain yield, Table S1: List of accessions, Table S2: Significant GWAS results for root related traits, Table S3: Significant GWAS results for grain yield, Table S4: Gene models identified within the 15 root QTL hotspots.

**Author Contributions:** Conceptualization, S.S. and J.M.S.; Data curation, R.R.; Formal analysis, R.R.; Funding acquisition, C.R. and J.M.S.; Investigation, R.R., S.S. and J.M.S.; Methodology, R.R., S.S. and J.M.S.; Project administration, C.R. and J.M.S.; Supervision, S.S. and J.M.S.; Writing—original draft, R.R.; Writing—review & editing, S.S., C.R. and J.M.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was funded by projects AGL2015-65351-R and RTA2015-00072-C03-01 of the Spanish Ministry of Economy and Competitiveness. R.R. is a recipient of a PhD grant from the Spanish Ministry of Economy and Competitiveness.

**Acknowledgments:** The authors acknowledge the contribution of the CERCA programme (Generalitat de Catalunya). Thanks are given to the group of Agricultural Genetics of University of Bologna for technical support.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

#### **References**

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## **Improvement of a RD6 Rice Variety for Blast Resistance and Salt Tolerance through Marker-Assisted Backcrossing**

#### **Korachan Thanasilungura 1, Sukanya Kranto 2, Tidarat Monkham 1, Sompong Chankaew 1 and Jirawat Sanitchon 1,\***

Received: 24 July 2020; Accepted: 30 July 2020; Published: 1 August 2020

**Abstract:** RD6 is one of the most favorable glutinous rice varieties consumed throughout the north and northeast of Thailand because of its aroma and softness. However, blast disease and salt stress cause decreases in both yield quantity and quality during cultivation. Here, gene pyramiding via marker-assisted backcrossing (MAB) using combined blast resistance QTLs (*qBl 1, 2, 11*, and *12*) and *Saltol* QTL was employed in solving the problem. To pursue our goal, the RD6 introgression line (RGD07005-12-165-1), containing four blast-resistant QTLs, were crossed with the Pokkali salt tolerant variety. Blast resistance evaluation was thoroughly carried out in the fields, from BC2F2:3 to BC4F4, using the upland short-row and natural field infection methods. Additionally, salt tolerance was validated in both greenhouse and field conditions. We found that the RD6 "BC4F4 132-12-61" resulting from our breeding programme successfully resisted blast disease and tolerated salt stress, while it maintained the desirable agronomic traits of the original RD6 variety. This finding may provide a new improved rice variety to overcome blast disease and salt stress in Northeast Thailand.

**Keywords:** gene pyramiding; aroma; QTL; chromosome; selection; introgression line

#### **1. Introduction**

Rice (*Oryza sativa* L.) is consumed as a staple food in Asia, especially in the southeast region. In Thailand, the indica rice variety RD6 developed from KDML105 through gamma irradiation is one of the most favorable glutinous rice consumed throughout the northeast of Thailand [1,2]. Because of its cooking quality, aroma, and softness, production demand has increased over time. However, its yield of 4.16 ton/ha fails to meet its potential, due to biotic and abiotic stress.

Rice blast disease caused by the fungus *Pyricularia grisea* (Cooke) Sacc. leads to crop losses up to 85% of total yield [3]. Disease symptoms occur in all stages of plant growth, beginning with blast discoloration and wilting of the foliage [4]. Neck blast can be found at the flowering stage, accelerating plant death [5]. Severe damage was also observed within areas of intensive planting with high doses of nitrogen application [6]. Development of new rice varieties resistant to blast fungus is an alternative approach to diminish or control the invasion of this pathogen. The resistance quantitative trait loci (QTL) have been investigated to achieve parental varieties, which are further used for gene pyramiding in breeding programmes. Currently, more than 100 blast-resistant genes have been identified, of which 22 genes structures have been cloned [7]. In Thailand, few studies of blast resistant genes have been conducted. Noenplab et al. [8] studied the relationship of leaf blast and neck blast of resistant genes in the Jao Hom Nil (JHN) variety, in which the resistant QTLs were detected on chromosomes 1 and 11. The resistant QTLs conferred resistance to both leaf blast and neck blast. Suwannual et al. [9]

pyramided four blast-resistant QTLs, individually, on chromosomes 2 and 12 within the P0489 variety, and on chromosomes 1 and 11 carried by the JHN variety, resulting in the creation of new RD6 introgression lines. Their results demonstrated that the RD6 introgression lines carrying a high number of QTLs (achieved through pyramiding) reached a broader spectrum of blast resistance to the blast pathogens prevalent in the region.

In addition to rice blast fungus, salt stress is a crucial constraint for RD6 production. Thailand's northeast region is the country's largest rice-producing area, and it is comprised of two basins: Sakon Nakhon and Nakhon Ratchasima. In those basins containing an understructure of accumulated salt rock, the salt-affected range covers approximately 1.84 Mha [10]. Evaporation during the dry season tends to raise salinity from the subsoil to the surface, thereby increasing salinity intensity and increasing salt stress from 2–4 dS/m to 8–16 dS/m [10–12]. Rice is a salt-sensitive crop, capable of tolerating salinity at moderate levels of electrical conductivity (4–8 dS/m) [13]. Therefore, rice produced under rain-fed, lowland conditions is usually exposed to high levels of soil salinity. The RD6 variety, which is well known, was identified as a geographical indication (GI) within the Tung Gula Rong Hai of the Northeast, Thailand. Specifically, RD6 requires optimal soil salinity to enhance rice seed aroma [14]. However, an abundance of salinity can reduce rice plant growth, tiller number, and seed set-up [15], and the stress caused by excessive salt can significantly reduce total crop yield and result in plant death [16].

The Pokkali variety, derived from the International Rice Research Institute (IRRI), has become a well-known source of salinity tolerance worldwide, attributed to the salt-tolerant QTL located on rice chromosome 1 (*Saltol*) [17–20]. Therefore, several researchers have attempted to develop salt-tolerant rice varieties using the *Saltol* QTL [21–26].

The marker-assisted backcrossing (MAB) method has been employed to obtain beneficial QTLs from donor parents via introgression between the qualitative and quantitative traits from landraces and wild relatives [27] due to the precision method with shortened time frame in both foreground and background selection. MAB provides effective gene selection and/ or QTLs for pyramiding multi-genes/QTLs within the rice population. These benefits further support breeding practices for improved resistance and tolerance [28–32]. The objective of this study was to determine the blast resistance and salt tolerance levels within the RD6 introgression lines by pyramiding four blast-resistant and one salt-tolerant QTL into the RD6 rice variety in both greenhouse and field conditions.

#### **2. Materials and Methods**

#### *2.1. Plant Materials and Marker-Assisted Backcrossing Selection (MABS)*

Three parental varieties/lines were used to generate the BC4F4 population, representing a pseudo-backcrossing approach to increasing the recurrent genetic background of a pyramiding population, comprised of the RD6 (recurrent parent), Pokkali (obtained for the *saltol* QTL present on chromosome 1), and RGD07005-12-165-1 (the RD6 near-isogenic line obtained from the Rice Gene Discovery Unit, Kasetsart University, Thailand). The RGD07005-12-165-1 obtained four blast-resistant QTLs from the JHN and P0489 varieties on chromosomes 1 and 11, and chromosomes 2 and 12, respectively. The breeding program was subsequently improved within the population through MAB, by crossing the RGD07005-12-165-1 with the Pokkali variety to improve salt tolerance. The F1 was then backcrossed with RGD07005-12-165-1through BC1F1, whereas BC1F2 was utilized as a marker-assisted selection (MAS) in the blast-resistant and salt-tolerant QTLs. In this step, the flanking marker RM3412/RM10748 was used for the selected *Saltol* QTL [33], RM319/RM212 and RM114/RM224 were used for selected blast-resistant QTLs on chromosomes 1 and 11, respectively [8], and RM48/RM207 and RM313/RM277 were used for selected blast-resistant QTLs on chromosomes 2 and 12, respectively [9], as shown in Figure 1.

**Figure 1.** Breeding schematics for the development and validation of the RD6 NILs populations.

Total genomic DNA from young leaves of individual plants, lines, and their parents was extracted according to the method described by Dellaporta et al. [34] with slight modifications. The PCR reactions for SSR markers were carried out in a volume of 10 μL, containing 25 ng of genomic DNA, 1 × PCR buffer, 1.8 mM MgCl2, 0.2 mM dNTP, 0.2 μM forward and reverse primer, and 0.05 unit Taq DNA polymerase. DNA amplification was performed in a DNA thermal cycle for five minutes at 95 ◦C, followed by 35 cycles of 30 s at 95 ◦C, 30 s at 55 ◦C, and two min at 72 ◦C, with a final extension of seven minutes at 72 ◦C. The amplification products were separated via 4.5% polyacrylamide gel electrophoresis [32]. The selected lines in the BC1, BC2, and BC3 generations were backcrossed with the RD6 variety, using MAS for *saltol* and blast-resistant QTL selection. Trait qualities; such as glutinous type, aromatic character, and gelatinization temperature (GT) in the BC2F2 populations were fixed through MAS using glutinous 23 primer on chromosome 6, badh2 on chromosome 8, and RM190 on chromosome 6 (Table S1). Each backcross generation within the BC4F4 populations was evaluated for salt tolerance and blast resistance (Figure 1).

#### *2.2. Evaluation of Salt Tolerance and Blast Resistance in the BC2F2:3 Populations (Exp. 1)*

The evaluation of salt tolerance and blast resistance of the BC2F2:3 lines, as well as the parental and check varieties, were conducted at the Department of Agronomy, Faculty of Agriculture, Khon Kaen University, Khon Kaen, Thailand. The salt tolerance evaluation was performed through two methods, salt solution and artificial soil salinity. The salt solution method was laid out in a completely randomized design (CRD), with four replications. Seedlings were transplanted at seven days to 50 × 57 cm2 Styrofoam sheets with 1.5 cm diameter holes [17]. Fertilizer was applied with

Yoshida nutrient solution [35] at three to twenty-one days of age. NaCl was then added to reach the electrical conductivity (EC) of 4 dS/m. EC was subsequently increased at three-day intervals, reaching EC 8 and 12 dS/m. When the susceptible checks (IR29) presented salt injury symptoms, salt tolerance data were recorded following the standard evaluation score (SES) [36]. The artificial soil salinity method was also laid out in CRD with four replications. Seedlings were transplanted at seven days and transferred from the spent soil trays to the water trays. At 14 days after transplanting, the experimental trays were treated with NaCl, which adjusted the water solution to EC 8dS/m, which increased to EC 12dS/m after two days. Salt tolerance data were recorded similarly in both methods.

Blast resistance was also evaluated via the upland short-row method at the Sakon Nakhon Rice Research Center, Sakon Nakhon, Thailand. The experiment was laid out in CRD with three replications. Seeds of each BC2F2:3 line were sown in rows (approximately) 50 cm long and 10 cm apart. A susceptible KDML105 variety was planted alternately with every two testing varieties. Blast resistance scores were recorded following the SES method [36].

#### *2.3. Evaluation of Salt Tolerance in the BC3F4 Populations (Exp. 2)*

The BC3F4 lines and parental varieties were evaluated for salt tolerance in field conditions. The experiment was conducted at the Ban Daeng Village, Ban Fhang, Khon Kaen, Thailand. The experiment was laid out in a randomized complete block design (RCBD) with three replications. Germinated seeds were sown on seedbeds; then, at thirty days, the seedlings were transplanted to the field. Plot sizes were 1 × 1.5 m2, in three rows, spaced 25 × 25 cm between and within rows. The RD6 variety was planted between every five plots within the test lines to ensure that salinity occurred uniformly. Fertilizer (23.44 kg/ha of N, P2O5, and K2O) was applied at four days after transplanting, and hand weeding and chemical application for disease and insect control were performed as needed. When the susceptible check (RD6) presented salt injury symptoms, salt tolerance data were recorded following SES [37]. Moreover, the agronomic traits including 1000/seed weight, seed length, seed width, and seed shape were recorded.

#### *2.4. Evaluation of Salt Tolerance and Blast Resistance Evaluations in the BC4F3 Populations (Exp. 3)*

The salt tolerance evaluation of the BC4F3 lines was conducted in greenhouse conditions at the Department of Agronomy, Faculty of Agriculture, Khon Kaen University, Khon Kaen, Thailand. The experiment was laid out in CRD with three replications, as described in the salt solution method of Exp. 1. The experiment for blast resistance of the BC4F3 lines was conducted at the Department of Agronomy, Faculty of Agriculture, Khon Kaen University, Khon Kaen, Thailand through field and upland short-row experiments. The upland short-row method was conducted similarly to Exp. 1. The field experiment was laid out in RCBD with three replications, in plots 0.5 × 1 m2, in three rows, spaced 25 × 25 cm between and within rows. Seedlings were transplanted at 30 days of age. Symptoms of natural blast infection were identified when seedlings began to show signs of infection, following the protocol of the SES [37], in which both leaf and neck blast symptoms were recorded.

Agronomic trait data, such as plant height (PH) and panicle length (PL), were recorded at pre-harvest; whereas post-harvest data recorded 4/panicle seed weight (SW4P), 1000/seed weight (1000SW), total dry weight (TDW), total seed weight (TSW), harvest index (HI), seed length (SL), seed width (SW), and seed shape (SS) (ratio of SL/SW). Additionally, seed qualities of the BC4F4 seeds, such as seed morphology, including SL, SW, SS and seed color of brown and paddy rice and aromatic traits, were evaluated and compared with the RD6 variety. The seed aromatic evaluation of each line was achieved through the quantitative determination of 2-acetyl-1-pyrroline (2AP) content using automated headspace gas chromatography following the methods as prescribed by Sriseadka et al. [38]. In brief, polished seed (1.00 g) were ground and then placed in a 20 mL headspace vial. The headspace vials were immediately sealed with PTFE/silicone septa and aluminum caps prior to analysis by static headspace-gas chromatography. A static headspace (Model 7697A, Agilent Technologies, Santa Clara, CA, USA) coupled to an Agilent 7890B Series GC system equipped with an Agilent

5977B GC/MSD system was used. A series of 2AP standard solutions with concentrations of 1.22, 2.45, 4.90, 9.79, and 19.67 ppm in isopropanol were prepared, which was added to headspace vials containing 1.00 g of non-aromatic rice seed (cv. Chai Nat 1) used as the external standard. The optimum headspace operating conditions were oven temperature 110 ◦C, loop temperature 120 ◦C, transfer line temperature 130 ◦C, vial equilibration time 10 min with high speed shaking, pressurizing time 0.15 min, loop equilibration time 0.40 min, and inject time 0.50 min. The headspace volatiles were separated using an HP-5 (25 m × 250 μm × 0.25 μm film thickness) column (J&W Scientific, Folsom, CA, USA). The optimum GC conditions were achieved using an HP-5 column with a splitless injection at 210 ◦C. The column temperature programme began at 50 ◦C and increased to 200 ◦C at 10 ◦C/min. Purified helium was used as the GC carrier gas at a flow rate of 1.2 mL/min. A calibration curve for 2AP analysis by headspace was generated by spiking known concentrations of 2AP into a non-fragrant rice variety (Chai Nat 1). Samples were run in triplicate, and the concentration of 2AP was calculated based upon the relative peak area of external standard.

#### *2.5. Evaluation of Salt Tolerance in the BC4F4 Populations (Exp. 4)*

The salt tolerance evaluations of the BC4F4 lines were conducted in greenhouse conditions at the Department of Agronomy, Faculty of Agriculture, Khon Kaen University, Khon Kaen, Thailand, via the salt solution method. Laid out in CRD with four replications, the planting methods and experiment protocol were similar to those described in Exp1, except that the EC in Exp. 4 was adjusted to 18 dS/m. The leaf, stem, root, and total dry weight data were recorded. Dry weight of the seedlings was determined by oven-drying the seedlings at 80 ◦C for 3 days, and the percentages of Na+ and K+ in the leaves, stems, and roots also were determined, according to the flame photometric method [39]. In addition, the percentages of Na+ and K+ in the leaves, stems, and roots also were determined, according to the flame photometric method [39]. In brief, pounding sample volume 0.5 g of each tissues were digest by 10 mL nitric acid (HNO3) and 5 mL perchloric acid (HClO4), then incubated at 200 ◦C. The contents were covered to reflux acid fumes generated during digestion until digest appeared translucent. After cooling down, 100 mL deionized water was added to each digestion tube. Contents were then vortexed and passed through qualitative cellulose filter paper (Whatman No.1, Sigma-Aldrich®, St. Louis, MO, USA) and measurement the K+ (768 nm wavelength) and Na+ (589 nm wavelength) by flame photometer (Model 410 Flame Photometer, Sherwood Scientific Limited, Cambridge, UK). The K+ and Na+ concentrations of samples were compared with the known standard solutions of 0.0, 5.0, 10.0, 15, and 20 ppm from a calibration curve with a correlation coefficient (r2) = 0.999. Finally, the amount of K+ and Na+ were transformed to percentage when compared with dry weight of raw materials.

#### *2.6. Data Analysis*

Salt tolerance scores, blast resistance scores, and agronomic trait data were analyzed via the STATISTIC 10© program (1985–2013) (Analytical Software, Tallahassee, FL, USA). Means were compared by the least significant difference (LSD) at *p* < 0.05.

#### **3. Results**

#### *3.1. Development of Populations Through MABS*

The F1 population (RGD07005-12-165-1 × Pokkali) was backcrossed with RGD07005-12-165-1 using MAS through BC2F2, in which the MAS contained genes of the glutinous type, aromatic, and gelatinization temperature (GT); and the BC2F2:3 populations were then evaluated for salt tolerance and blast resistance. One BC2F2:3 line (no. 74) obtained all QTLs required to obtain the target traits necessary for generation advancement (Table S2). In developing the BC3 generation, BC2F2:3 (no. 74) was crossed back to RD6 to produce the BC3F1, in which the BC3F1 was developed through the BC3F4 population by MAS. Thirty-one BC3F4 lines were evaluated for salt tolerance and agronomic traits

within the salted field. Among the BC3F4 lines, the BC3F4 (no. 132) demonstrated salt tolerance and agronomic characteristics similar to those of the RD6 variety and were subsequently selected for the development of the BC4F1 through BC4F3. The BC4F3-132-12-61 line was successfully developed via the MABS, consisting of one QTL for salt tolerance and four QTLs for blast resistance. BC4F3 lines with varied QTL combinations were also validated in the tests.

#### *3.2. Evaluation of Salt Tolerance and Blast Resistance in the BC2F2:3 Populations*

Eight of the BC2F2:3 lines were evaluated for salt tolerance in the seedling stage through the salt solution and artificial soil salinity methods. The results concluded that the BC2F2:3 lines presented as highly significant at EC 12 dS/m. BC2F2:3 (nos. 23, 67, and 74), demonstrated tolerance (T) in both methods, similar to that of the tolerant check (Pokkali); while the remaining five BC2F2:3 lines showed moderate tolerance (MT) under both the salt solution and artificial soil salinity tests (Table 1). The BC2F2:3 lines were evaluated for blast resistance via the upland short-row method, through which the blast resistance reaction (R) was similar to that of the resistance checks (JHN and P0489). Their resistance proved greater than that of the susceptible check (RD6), due to the presence of blast-resistant QTLs, and superior to that of the original RD6 (Table 1).

**Table 1.** Salt tolerance scores at EC 12 Ds/m in salt solution, artificial soil salinity, and blast resistance within the BC2F2:3 populations.

T = tolerance, MT = moderate tolerance, S = susceptible, HS = highly susceptible, R = resistant, MR = moderate resistance, MS = moderately susceptible. \*\* = significant different at *p* < 0.01. Different letters after the mean within a column showed a significant difference. CV = the coefficient of variation.

#### *3.3. Evaluation of Salt Tolerance in the BC3F4 Populations*

BC2F2:3 (No.74) was selected for the next backcross cycle with the RD6, and MAS was performed through the BC3F4 population. Thirty-one BC3F4 lines were screened for salt tolerance and agronomic performance. The results showed that the 21 BC3F4 lines were salt-tolerant (T), whereas the 10 BC3F4 lines proved only moderately tolerant (MT). Agronomic traits of the BC3F4 lines were similar to those of the RD6 variety, which was in accordance with our objectives (Table 2). Five BC3F4 lines (nos.22, 36, 115, 129, and 132) demonstrated superior tolerance (T) and agronomic traits (1000/SW, SL, SW, and SS), again, similar to those of the original RD6, and were selected as donor parents for the development of the BC4. The results showed that BC3F4 (no. 132) produced the best performance for pollination, and was therefore selected for development of the BC4 population.

STS = salt tolerance score, 1000SW = 1000 seeds weight, SL = seed length, SW = seed width, SS = seed shape (ratio of SL/SW). = significantly different at *p* < 0.01. Different letters the mean within a column showed a significant difference. CV=the coeffcient of variation.

#### *3.4. Evaluation of Salt Tolerance and Blast Resistance in the BC4F3 Populations*

This experiment evaluated eight BC4F3 lines obtained via *Saltol* QTL from MAS, specifically, the combination lines, donor parents, and tolerance, and susceptible checks were screened for salt tolerance via the salt solution method (EC 12 dS/m). The results showed highly significant scores within the BC4F3 population and the tolerant check (Pokkali), with salt tolerance scores less than 5.0, indicated as tolerant (T). The recurrent parent (RD6) presented moderate tolerance (MT) with a score of 6.5 (Table 3). The results indicated that the Saltol QTL within the BC4F3 population demonstrated superior performance in salt tolerance to that of the recurrent parent (RD6).

The BC4F3 populations were evaluated for blast resistance in field conditions through natural infection. Since blast infection in the field was not severe, leaves with severe symptoms from the surrounding trap plants were collected and stored in a bag, under dark conditions, for twelve hours to induce spore formation. The natural inoculum was added with water, and spore suspension was then sprayed over the field. The BC4F3 populations and their parents showed a high level of blast resistance. A total of eight BC4F3 lines demonstrated high resistance (HR), similar to those of both the donor parents, whereas the recurrent parent (RD6) presented moderate susceptibility (MS) (Table 3). Importantly, some introgression lines were greater in resistance than the recurrent parent (RD6). The second peak of bimodal rain, which occurred in the flowering stage, initiated significant signs of neck blast. The results found that eight BC4F3 lines showed resistance (R) similar to that of the donor parent, whereas the RD6 was moderately susceptibility (MS) to neck blast (Table 3).

Additionally, blast resistance evaluation was also conducted via the upland short-row method in the seven BC4F3 lines, which presented highly significant blast-resistant levels. The BC4F3 lines showed resistance abilities similar to both donor parents, whereas the RD6 recurrent parent presented moderate resistance (MR). Notably, the blast-resistant genes in the BC4F3 populations provided blast resistance in the seedling stage similar to that of the tilling and grain filling stage, and evidenced greater resistance than that of the RD6 recurrent parent (Table 3).

The agronomic traits were also evaluated in the blast field experiments. The results were highly significant within the BC4F3 populations and their parents for ten traits: PH; PL; SW4P; 1000/SW; TDW; TSW; HI; SL; SW; and SS (Table 4). The BC4F3 132-12 maintained agronomic traits similar to those of the recurrent parent (RD6) for nine traits, except for the 1000/SW (1000 seed weight), in which the RD6 presented moderate susceptibility (MS) for leaf and neck blast, resulting in low grain filling (Table 4). The results indicate that the newly developed MAB population had greater resistance than that of the original RD6 variety while maintaining its desirable agronomic traits and satisfying consumer demand. The BC4F3 lines also showed improvements in seed length (SL), seed width (SW), and seed shape (SS) as long and slender type. Additionally, the color of paddy rice was straw yellow, similar to the original RD6 varieties (Figure 2).

HR = high resistance, R = resistance, MR = moderate resistance, MS = moderate susceptible, S = susceptible. \*\* = significantly different at *p* < 0.01. Different letters after the mean withina column show a significant difference. CV=the coefficient of variation.

#### *Agronomy* **2020** , *10*, 1118

column show a significant difference. CV = the coefficient of variation.

**Figure 2.** Seed quality, seed length, seed shape, and seed color of 10 seeds of the BC4F3 populations compared with the RD6 variety. (**a**) RD6, (**b**) BC4F3 133-12-61, (**c**) BC4F3 133-25-18, (**d**) BC4F3 132-14, (**e**) BC4F3 132-98, (**f**) FC4F3 132-174, (**g**) BC4F3 132-51, (**h**) BC4F3 132-167, (**i**) BC4F3 132-276. The top row is brown rice and bottom row is paddy rice with straw color. The small scale is in millimeters.

Eight BC4F3 lines and the RD6 (recurrent parent) were evaluated for seed aroma through the determination of 2AP content via the automated headspace gas chromatography method. The results showed a significant difference among the BC4F3 lines and RD6 variety, with mean values of 2AP content of the BC4F3 exceeding 3.00 ppm (Table 4). BC4F3 132-98-87 presented the highest 2AP content (4.68 ppm), similar to that of the RD6 (Table 4). The BC4F3 lines and RD6 were determined to be similar in fragrance.

#### *3.5. Evaluation of Salt Tolerance in the BC4F4 Population*

Two BC4F4 lines, BC4F4 132-12 and BC4F4 132-167, presented as tolerant (T) in salt tolerance evaluations, similar to that of the Pokkali, whereas the remaining five BC4F4 lines and the recurrent parent (RD6) showed moderate tolerance (MT) (Table 5). Statistically, seven of the BC4F4 lines demonstrated EC values up to 18 dS/m. The dry weights of the seven BC4F4 lines, their parent, and KDML105 check varieties presented as highly significant, in which the tolerant check (Pokkali) presented the highest LDW, SDW, RDW, and TDW evident in the weights of leaf stems and roots, whereas the BC4F4 lines were similar to that of the recurrent parent (RD6) (Table 5). The results indicate that salinity significantly affected dry weight.

**Table 5.** Leaf, stem, root, and total dry weights in the experiment 4.

LDW = leaf dry weight, SDW = stem dry weight, RDW = root dry weight, TDW = total dry weight. \*\* = significantly different at *p* < 0.01. Different letters after the mean within a column show a significant difference. CV = the coefficient of variation.

Additionally, the percentages of Na+ and K+ in leaf stems and roots were also recorded and proved highly significant in all traits. The tolerant checks (Pokkali) presented the lowest Na+ in leaves (1.94), stems (1.75), and roots (4.94), whereas K+ in leaf stems were non-significant within the BC4F4 lines. Moreover, Pokkali also displayed low Na+-to-K+ ratios in rice shoots. The BC4F4 132-12 also showed low levels of Na+ in both stems and leaves, as well as low Na+-to-K+ ratios in rice shoots comparable to those of the Pokkali variety, yet lower than those of the RD6. We may infer that transferring the saltol QTL from the Pokkali line created greater tolerance in subsequent breeding lines that that of RD6 recurrent parent (Figure 3). Interestingly, the BC4F4 lines presented Na+ levels in leaves and stems similar to those in Pokkali, and significantly different in RD6. This confirmed the Pokkali saltol QTL's ability to exclude Na+ from the leaf blade, expressing salt tolerance through the salt-excluder method. Moreover, the salt tolerance scores were negatively correlated to the leaf, stem, root, and total dry weights (r = −0.7899 \*\*, r = −0.8136 \*\*, r = −0.8140 \*\*, r = −0.8065 \*\*, respectively) and positively correlated with leaf and stem Na+ (r = 0.7670 \*\* and r = 0.8917 \*\*, respectively) (Table 6). The results indicated that salt stress decreased plant growth in susceptible varieties, while growth was maintained in salt-tolerant varieties and that the accumulation of Na+ in leaves and stems was related to salt susceptibility (Table 6).

**Table 6.** The correlation between leaf dry weight, stem dry weight, root dry weight, total dry weight, leaf Na+, leaf K+, stem Na+, stem K+, root Na+, root K+, and the salt score of seeds in the BC4F4 populations and check varieties in the experiment 4.

LDW = leaf dry weight, SDW = stem dry weight, RDW = root dry weight, TDW = total dry weight. \* = significantly different at *p* < 0.05. \*\* = significantly different at *p* < 0.01.

**Figure 3.** Na+ and K+ on leaves, stems, and roots in the experiment 4. Different letters within a color bar show a mean significant difference of each line.

#### **4. Discussion**

Since its release in 1977, the RD6 glutinous rice variety has remained a staple food crop for domestic consumption in Thailand's north and northeast regions. Comprising 83% of total glutinous rice production in these areas, consumers have developed a preference for its superior characteristics. However, the RD6 variety suffers from several production constraints, including biotic stress responsible for both rice blast [40] and bacterial blight disease [41]. Current research has attempted to eliminate sustainable infection-resistant production practices by pyramiding multiple resistant genes [42]. To date, an RD6 introgression line capable of resisting both biotic and abiotic stress has yet to be developed. Thailand's salt rock basins of Sakon Nakhon and Nakhon Ratchasima have demonstrated that consistent levels of salinity can enhance the fragrance of the RD6 rice variety [14] and increase production.

This study proposes the successful introgression of blast-resistant QTLs (*qBL 1, 2, 11*, and *12*) from RGD07005-12-165-1 and the *Saltol* QTL (Pokkali, chromosome 1) to improve the RD6 rice variety through the MAB method within BC4F4 populations. Trait evaluations were completed for the validation of progenies with desirable traits in each advanced population, based on the introgression of the genetic foregrounds and maintenance of the genetic backgrounds, respectively.

Salt salinity was absent in several areas of northeast Thailand, due to high levels of NaCl [43], and such factors as precipitation, soil type, and field management. In past research, the evaluation of salt tolerance was typically conducted through salt screening, hydroponic culture, and soil culture, as well as through pot and field methods [33]. The current study assessed salt tolerance within the breeding populations studied through salt solution, artificial salt culture, and field condition evaluations. Based on the results, salt evaluation under field conditions produced the lowest capability among the tested rice lines (Tables 1–3 and 5), due to the inherent difficulties and uncertainties present under field conditions. Kranto et al. [33] reported that effective alternative screening approaches must be proven to correlate with results produced within the early phases of growth in both greenhouse and field conditions. Within the present study, visual symptom scores of salt stress generated through the salt solution method proved to be the most appropriate method with which to confirm tolerance abilities within a breeding population (Table 5), suggesting that the salt solution method could, therefore, substitute salt tolerance score analysis in field conditions (Table 5). However, we acknowledge the necessity to evaluate RD6 plant types, yield performances, and agronomic traits within the field. The introgression lines developed within our study were evaluated for similarity with the original RD6 agronomic traits, namely plant height, panicle length, 4/panicle seed weight, 1000/seed weight, total dry weight, total seed weight, harvest index, seed length, seed width, and seed shape, as well as seed qualities, such as seed morphology (Tables 2 and 4).

The *Saltol* QTL on chromosome 1 from the Pokkali rice variety has been commonly used for rice improvement in several studies [21–26]. In our results, the *Saltol* QTLs from the Pokkali variety produced the greatest salt tolerance within the RD6 introgression lines (Table 1, Table 2, Table 3, and Table 5). This *saltol* QTL also contributed to the maintenance of low Na+, high K+, and low Na+/K+ homeostasis levels in rice stems, further resulting in increased salt tolerance [24,44] (Figure 3). The Pokkali variety was classified to balance the influx of Na+ and K+ for dilution in the mechanism, creating the ability to exclude Na+ from leaf blades and stems [45,46]. As the water up-take mechanisms in rice accept both nutrients and salt together, the Pokkali variety thereby demonstrated the highest and most significant differences in leaf, stem, root, and total dry weights when compared with other breeding lines (Table 5). However, the BC4F4 lines presented the agronomic traits (above) more closely matched to the RD6 than to the Pokkali (Table 5), due to the advance generation and visual selection of the trait performances (Tables 2 and 4). RD6 performance is very important for farmer acceptance and crop adaptation in our test areas. For example, excessively tall RD6 rice plants present problems in the grain filling stages as a result of heavy wind or rain [47]. Visual selection may explain the differences in percentages of Na+ of the RD6 introgression lines with those of Pokkali (Table 5, Figure 3).

As a photosensitive rice variety, the RD6 grows once a year, during Thailand's rainy season from late May to November [48]. These bimodal rain patterns produce favorable conditions for the occurrence of blast disease, causing damage in all stages of growth. Leaf blast generally occurs during the seedling and tilling stages, whereas neck blast usually occurs during the reproductive phase [4]. In our study, introgression lines were evaluated for blast disease in both the field and upland short-row evaluations.

In this study, the upland short-row method displayed greater incidences of blast disease, due to the favorable microclimate and moisture contents around the experimental plots (Tables 1 and 3) [49]. The experimental field was influenced by bimodal rain, capable of inducing leaf and neck blast symptoms (Table 3), further indicating the resistance of the QTLs [42,50]. Noenplab et al. [8] also reported that the blast QTL on chromosome 11 in the JHN variety successfully contributed to leaf and neck blast resistance. Pyramiding of four blast-resistant QTLs through MAS achieved high levels of blast resistance and broad-spectrum resistance to pathogens prevalent in the region [9]. Moreover, the testing of RD6 introgression lines for durable blast resistance and no-yield penalties were observed [42]. The results, herein, further demonstrated that neck blast disease caused direct yield loss during the grain filling phase [51], as well as lower 1000/SW within the original RD6 variety compared with those of the RD6 introgression lines (Table 4).

The resistance/tolerance abilities of the RD6 introgression lines represent the foreground genetics capable of enhancing plant breeding programs. However, maintaining the background of the original RD6 variety is also desirable; therefore, the quality and performance of the RD6 within the QTL introgression was also a consideration. The BC4F4 populations, herein, were achieved through the introgression of blast-resistant QTLs (*qBL 1, 2, 11*, and *12*) from RGD07005-12-165-1 and *Saltol* QTL (Pokkali) and improved the RD6 rice variety through MAB. Consequently, the performance of the RD6 introgression lines was similar to that of the original RD6 variety (Table 4, Figure 2). The results indicate that foreground and background selection, together with visual selection, accurately depicts the efficiency of MAB.

#### **5. Conclusions**

Improvement of the RD6 rice variety for salt tolerance and blast resistance was successfully achieved utilizing the *Saltol* QTL and *qBl* (*1*, *2*, *11*, and *12*) through marker-assisted backcrossing, together with phenotypic selection. The resulting BC4F4 132-12 introgression line exhibited superior salt tolerance, blast resistance, and reduced neck blast and was capable of maintaining higher qualities and agronomic performances than that of the original RD6 variety.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4395/10/8/1118/s1, Table S1: QTL traits and primer sequences of the SSR markers for blast resistance and salt tolerance, Table S2: Genotype of the BC2F2 populations derived from MAB.

**Author Contributions:** K.T. and S.K. conceived the study. T.M., S.C., and J.S. designed the experiments. K.T. and S.K. performed the experiments. J.S. and S.C. supervised the study. K.T., T.M., S.C., and J.S. wrote the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** This research was supported by The Plant Breeding Research Center for Sustainable Agriculture and The Salt-Tolerance Rice Research Group, Khon Kaen University, Khon Kaen, Thailand. Our gratitude is also extended to the Sakon Nakhon Rice Research Center for their support in our field experiments and Ubon Ratchathani Rice Research Center for their support in 2AP analysis.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Identification of QTLs Controlling Resistance**/**Tolerance to** *Striga hermonthica* **in an Extra-Early Maturing Yellow Maize Population**

#### **Ba**ff**our Badu-Apraku \*, Samuel Adewale, Agre Patern, Melaku Gedil and Robert Asiedu**

International Institute of Tropical Agriculture (IITA), PMB 5320, Ibadan 200001, Nigeria;

<EMAIL_ADDRESS>(S.A<EMAIL_ADDRESS>(A.P<EMAIL_ADDRESS>(M.G<EMAIL_ADDRESS>(R.A.)

**\*** Correspondence<EMAIL_ADDRESS>Tel.: +234-8108482590

Received: 12 July 2020; Accepted: 6 August 2020; Published: 10 August 2020

**Abstract:** *Striga hermonthica* parasitism is a major constraint to maize production in sub-Saharan Africa with yield losses reaching 100% under severe infestation. The application of marker-assisted selection is highly promising for accelerating breeding for *Striga* resistance/tolerance in maize but requires the identification of quantitative trait loci (QTLs) linked to *Striga* resistance/tolerance traits. In the present study, 194 F2:3 families of TZEEI 79 × TZdEEI 11 were screened at two *Striga*-endemic locations in Nigeria, to identify QTLs associated with *S. hermonthica* resistance/tolerance and underlying putative candidate genes. A genetic map was constructed using 1139 filtered DArTseq markers distributed across the 10 maize chromosomes, covering 2016 cM, with mean genetic distance of 1.70 cM. Twelve minor and major QTLs were identified for four *Striga* resistance/tolerance adaptive traits, explaining 19.4%, 34.9%, 14.2% and 3.2% of observed phenotypic variation for grain yield, ears per plant, *Striga* damage and emerged *Striga* plants, respectively. The QTLs were found to be linked to candidate genes which may be associated with plant defense mechanisms in *S. hermonthica* infested environments. The results of this study provide insights into the genetic architecture of *S. hermonthica* resistance/tolerance indicator traits which could be employed for marker-assisted selection to accelerate efficient transfer host plant resistance genes to susceptible genotypes.

**Keywords:** maize (*Zea mays* L.); *Striga* resistance/tolerance; QTL mapping; F2:3 biparental mapping; Marker-assisted selection

#### **1. Introduction**

Maize (*Zea mays* L.) is the most widely grown staple food crop in sub-Saharan Africa (SSA), and accounts for a large proportion of carbohydrates, proteins, lipids and vitamins for millions of people in the sub-region [1,2]. The root hemi-parasitic plant *Striga hermonthica* is an important biotic constraint limiting maize production in SSA. The *S. hermonthica* problem in SSA is the result of the shift from traditional cereal based farming system which facilitated longer fallow periods which ensured that the soil *Striga* seed bank was maintained at levels that plants could tolerate [1]. However, the pressure on agricultural land has necessitated land use intensification, cereal mono-cropping and reduced fallow periods resulting in increased *Striga* seed bank and infestation levels that threaten the livelihood of millions of farmers [3,4]. De Groote et al. [5] reported that over six million hectares of agricultural land in Western, Eastern and Southern Africa are seriously affected by *Striga*. Reduction in grain yield due to *S. hermonthica* parasitism may be up to 100% under severe infestation and unfavorable environmental conditions such as low soil fertility, erratic rainfall patterns and low-input conditions [6–8]. The subsistence farmers are usually the most severely affected. The extent to which *S. hermonthica* affects the growth of its host varies tremendously, depending on the level of host plant resistance/tolerance, extent of infestation, and the prevailing environmental

conditions [9,10]. Resistance to *Striga* denotes the capability of the host plant to induce the germination of *Striga* seeds but prevents the parasite from attaching to the roots of the maize plants or kills the attached parasitic plants. Under *S. hermonthica* infestation, the resistant genotype supports considerably fewer *Striga* plants and produces a greater yield than the susceptible genotype [11–13]. Contrarily, a *Striga* tolerant genotype supports as many *Striga* plants as the sensitive or susceptible genotype [14] but produces more dry matter and shows fewer damage symptoms [15]. *Striga* damage in maize is used as the indicator of tolerance while emerged *Striga* plants is the indicator of resistance. The identification of maize genotypes that combine outstanding levels of resistance and tolerance is a promising breeding strategy and has been recommended for *Striga* resistance breeding in several studies [12–14,16,17]. In selecting for tolerance or resistance and high grain yield under *Striga* infestation, the primary traits of interest are the *Striga* damage and number of emerged *Striga* plants. Presently, maize genotypes with combined resistance and tolerance to *S. hermonthica* (possessing both low *Striga* damage and few *Striga* emergence counts) as well as high grain yield have been identified in the International Institute of Tropical Agriculture - Maize Improvement Program (IITA-MIP). *Striga* damage rating score is positively associated with *Striga* emergence counts, and the two traits are negatively associated with yield under *S. hermonthica* infested conditions. Similarly, large positive additive genetic correlation was recorded between grain yield and ears per plant as well as moderately large negative genetic correlations between grain yield and flowering traits [18]. Comparable results were reported by earlier workers [19,20]. Nevertheless, the genotypic correlation between *S. hermonthica* damage rating and *S. hermonthica* emergence counts have been found to be low, implying that different genes control the inheritance of the two traits [18,21].

*Striga* infestation is dependent on *Striga* seedbank in the soil resulting from continuous cropping of host plants, leading to the accumulation of *Striga* seeds which can remain dormant in the soil for more than a decade [22]. The germination of *Striga* seed is induced by the production of plant hormones called strigolactones produced by the maize plant in the roots. The hormones are released when the plant is under stress [23]. For the germinated *Striga* seed to survive as an obligate parasite, it must produce haustoria that attach to the roots of the maize plant through which it draws water and photosynthates [24]. Even though *Striga* possesses chlorophyll for photosynthesis, it still depends on its host for survival [25], as a result of its inability to accumulate enough photosynthates for autotrophic growth. Therefore, *Striga* establishes direct xylem links with the root system of the host [26] to obtain nutrients from its xylem sap [27]. Furthermore, a higher rate of transpiration in *Striga* compared to the host plant creates a water potential gradient from the host to the *Striga* plant [28]. *Striga* phytotoxic effects on the maize plant direct the partitioning of assimilates into the roots rather than the shoots for grain filling, thereby resulting in plant biomass and yield reduction [29].

Presently, management strategies of *S. hermonthica* include cultural, chemical, and biological approaches, which are non-economical and/or knowledge-intensive for subsistence farmers [6]. Planting of *Striga*-resistant maize varieties is presently considered the best control strategy and easy to adopt or deploy, particularly when combined with other management practices [30–33]. Resistance to *S. hermonthica* parasitism is mainly attributed to low production of *Striga* germination stimulants by the host plant, attachment of few *Striga* plants to the roots of the host plant as well as fewer *Striga* emergence [13,34]. When breeding for *S. hermonthica* resistance in maize, a combination of these resistance mechanisms is desirable in achieving effective and durable resistance [33]. The slow rate of development and deployment of *Striga* resistant genotypes is largely attributable to the complex genetics of resistance as well as limited knowledge of the specific mechanisms associated with resistance to *Striga* [35]. The resistance to *S. hermonthica* in maize is regulated by many genes or quantitative trait loci (QTL) with small additive effects and it is significantly influenced by the environment [13,36]. Therefore, breeding for *Striga*-resistant cultivars using conventional approaches by selecting maize cultivars with enhanced resistance which requires evaluation in multi-locations and years, has been less effective and time-consuming [27].

Marker-assisted breeding makes use of genotypic data in the identification of genotypes possessing desirable alleles, using linked genetic markers. Breeders employ marker-assisted selection (MAS) when an important trait that is difficult to assess phenotypically, is tightly linked to a molecular DNA marker that can be scored quickly and precisely [37]. QTL mapping approaches are important genomic tools employed in dissecting the genetic architecture of complex traits [38,39] as well as identification of genetic linkage through wide genotyping of a panel of germplasm displaying contrasting phenotypes across different environments [40]. For QTL identification, the development of next-generation sequencing technology has become a practicable technique to rapidly identify large number of single nucleotide polymorphisms (SNPs) throughout the genome [41]. Unlike the use of second-generation molecular markers which result in low-quality mapping, SNP markers provide new insights to rapidly identify QTL of interest.

Information on map positions of genes and linked markers on a chromosome are crucial for efficient determination of genetic architecture of polygenic traits in crop plants [42]. Several QTLs and candidate genes controlling resistance to *Striga* have been reported in cereals. For instance, Swarbrick et al. [32] identified three *S. hermonthica*-resistant QTLs in Kasalath–Koshihikari rice backcross inbreds, two of these QTLs originated from the Kasalath allele and one from the Koshihkari allele. The largest-effect QTL (Kasalath-derived allele) explained 16% phenotypic variance in the mapping population and was located on linkage group 4. Haussmann et al. [43] detected molecular markers associated with *S. hermonthica*-resistant QTLs, with the most significant QTL corresponding to the major-gene locus low germination stimulant (LGS) in linkage group I. Five genomic regions (QTLs) linked to stable *S. hermonthica*-resistant alleles from resistant variety N13 were detected through evaluation across a large number of field trials in Mali and Kenya. However, limited reports are available on the QTLs and genes controlling *Striga* resistance in maize. In a recent study by Adewale et al. [44] to identify molecular markers linked to *S. hermonthica* resistance in maize, 24 SNPs significantly associated with *S. hermonthica* resistance indicator traits under artificial *S. hermonthica* infestation were detected. The authors also identified four candidate genes on chromosomes 3, 5, 9 and 10, with functions closely associated with maize plant defense mechanisms against *S. hermonthica* parasitism. Identification of QTLs linked to *S. hermonthica* resistance followed by gene introgression into elite genetic backgrounds, has the potential to reduce yield losses due to *Striga* and will ultimately provide solid foundation for improving *Striga* resistance [45].

The objectives of this study were to identify QTL and underlying candidate genes conferring resistance/tolerance to *S. hermonthica* in maize using F2:3 biparental mapping population derived from a cross between a *Striga* resistant inbred line, TZEEI 79 and *Striga* susceptible inbred TZdEEI 11.

#### **2. Materials and Methods**

#### *2.1. Germplasm and Phenotyping*

Based on the reports of previous studies, two extra-early maturing yellow inbred lines, TZEEI 79 (*Striga* resistant/tolerant) and TZdEEI 11 (*Striga* susceptible) were selected as parents to generate the F2:3 progenies used in the present study [1]. TZEEI 79 is an outstanding *S. hermonthica* resistant/tolerant, drought and low-soil N tolerant inbred line developed in the IITA-MIP from the broad-based *S. hermonthica* resistant/tolerant as well as drought and low-soil N tolerant population, TZEE-Y Pop STR C0. TZEEI 79 has significant positive GCA (general combining ability) for grain yield as well as significant negative GCA effects for *Striga* damage and *Striga* emergence counts under *Striga* infestation and has been extensively used in the IITA hybrid program as an important resource for developing high-yielding, multiple stress tolerant hybrids as well as an efficient tester for classifying other inbreds into heterotic groups [46]. Crosses were made between TZEEI 79 and TZdEEI 11 designated as P1 and P2 respectively, to generate 220 F1 progenies. The F1 progenies and the parental lines were planted, and leaf samples were collected at 3 weeks after planting. Verification of the parental type alleles (quality control analysis) was carried out on the F1 progenies prior to advancement to F2. The F1 progenies were screened using two SSR primers (bnlg 182 and umc 1568) which were found to be polymorphic between the two parents. The analysis identified 170 true-to-type F1 hybrids which were advanced to F2. The 170 true-to-type F2 ears were planted ear-to-row and 194 F2 individuals which were randomly selected and selfed were used in the present study.

The F2:3 progenies and the two parental lines were screened under artificial *S. hermonthica* infestation at Mokwa (9◦18 N, 5◦4 E, 210 m above sea level, 1100 mm yearly rainfall, luvisol soil) and Abuja (9◦16 N, 7◦20 E, 445 m above sea level, 1500 mm yearly rainfall, ferric-luvisol soil) in the Southern Guinea savanna of Nigeria in 2018. At each experimental site, the trial was laid out using randomized incomplete block design (14 × 14 lattice) with two replicates. The experimental units were 3 m long single-row plots, with an inter-row spacing of 0.75 m and within-row spacing of 0.4 m, to achieve a target population density of 66,666 plants/ha. The fields for artificial *S. hermonthica* infestation at Mokwa and Abuja were treated with ethylene gas at 2 weeks before planting to eliminate any potential *Striga* seeds present in the soil. The *S. hermonthica* seeds used for the experiment were obtained from sorghum farms around the test locations at Abuja and Mokwa in 2017. The artificial *S. hermonthica* infestation was carried out as proposed by the IITA Maize Program [16]. Briefly, about a week before inoculation, the *S. hermonthica* seeds were carefully mixed with finely sieved sand at the ratio 1:99 by weight to ensure rapid and uniform infestation. A standard scoop calibrated to deliver approximately 5000 germinable seeds per hill was utilized for the artificial infestation. Three maize seeds were planted per infested hill and the seedlings were later thinned to two plants per stand at 2 weeks after emergence. Fertilizer application on the maize plots was delayed till about 30 days after planting, in order to subject the maize plants to stress, a condition that was expected to enhance strigolactone production. This ensured good germination of *Striga* seeds and attachment of *Striga* plants to the roots of host plants. At this plant growth stage, 20–30 kg Nha−1, 30 kg each of P and K were applied as NPK 15-15-15, taking into consideration the fertility status of the soil. Reduction in the fertilizer application rate was important because *Striga* emergence decreases at high N rate [16]. At 10 weeks after planting, typical symptoms of *Striga* infestation on the host plants were observed, such as chlorosis, leaf scorching (firing) and blotching, stunting, decrease in ear and tassel size, brown necrotic spots, leaf wilting and rolling, stalk lodging, open-tip of ears at late growing stage, and premature death of host plants. Host plant *Striga* damage severity was scored using a scale of 1 to 9. Rating scales 1–5 indicated resistance while 6–9 indicated susceptibility, where 1 = normal plant growth, no obvious symptoms, and 9 = all leaves completely scorched, collapse of host plants and no ear formation [16]. In addition, data were collected on *Striga* emergence count at 10 WAP as the number of *Striga* plants thriving on the maize root system as well as ears per plant (EPP) by dividing the total number of ears harvested per plot by the number of plants in a plot at harvest. Grain moisture was determined using Kett moisture tester PM-450 and grain yield (kg/ha) was calculated using the field weight of harvested ears per plot, adopting a shelling percentage of 80, adjusted to 15% moisture content [47].

#### *2.2. DArTseq Genotyping and SNP Data Filtering*

Young and healthy leaf samples from single plants of the F2 individuals and bulk samples from the parental lines were collected and frozen immediately after harvesting using liquid nitrogen and thereafter stored at −80 ◦C. Genomic DNA extraction was carried out following the DArT protocol (www.diversityarrays.com/files/DArT\_DNA\_isolation/). The extracted DNA was assessed for quality by visualization on agarose gel (2% w/v) and the quantity was estimated on NanoDrop-1000 spectrophotometer (NanoDrop, Wilmington, DE, USA) using the absorbance ratio A260/A280 to determine the concentration (ng/μL) and purity level of the DNA.

Genotyping of the 194 F2 individuals plus the two parents was carried out using DArTseq technology [48,49]. Genome complexity reduction which involved the use of a combination of two restriction enzymes (*Pst*I–*Mse*I) was used to create a genome representation of the analyzed samples. All fragments generated were amplified and sequenced to identify the single nucleotide polymorphisms (SNPs) using a proprietary analytical pipeline developed by DArT P/L. After a strict quality control process, which included parameters such as call rate, data reproducibility (~20% of samples replicated), and rate of monomorphism to eliminate monomorphic markers, 9951 SNPs were extracted from the evaluated germplasm. The 9951 SNPs were filtered for unmapped markers, duplicate markers and markers segregating between the *Striga* resistant and *Striga* susceptible parents. A total of 1139 high-quality DArTseq markers distributed across the 10 maize chromosomes were retained for the construction of genetic linkage map as well as the QTL mapping.

#### *2.3. Data Analysis*

The data recorded on emerged *Striga* counts and *Striga* damage severity scores were subjected to natural logarithm transformation. Thereafter, data collected on grain yield, ears per plant, *Striga* damage as well as *Striga* emergence counts were tested for normality using Shapiro–Wilk's (*W*) test [39,50] before analysis of variance. Box plots were made to visualize the distributions of grain yield and other traits under each research environment using ggplot2 library [51]. Analysis of variance was conducted across research environments using the general linear model procedure (PROC GLM) implemented in the Statistical Analytical System (SAS), version 9.3 [52]. In the analysis, environment, replications (environments), blocks (replications × environments) were considered as random and the F2:3 families (genotypes) as fixed effects. Estimates of broad sense heritability of the traits (H2) across research environments were computed on a family-mean basis as proposed by Holland et al. [53], using the following formula:

$$\mathcal{H}^2 = \frac{\sigma\_{\mathcal{S}}^2}{\sigma\_{\mathcal{S}}^2 + \frac{\sigma\_{\mathcal{S}^c}^2}{\varepsilon} + \frac{\sigma\_{\mathcal{S}}}{rc}}$$

where σ2 *g* = variance component due to the genotypes, σ2 *ge* = genotype × environment variance, σ*e* = experimental error variance; *e* = number of environments, and *r* = number of replications within environment.

Correlation coefficients were estimated among the traits with the adjusted means of the F2:3 families using the Ggally function implemented in GGally package [54]. Furthermore, the mixed linear model (MLM) established in META-R software [55] was used to compute the best linear unbiased estimates (BLUEs) for each genotype in each and across environments which were for the QTL analysis. The R/qtl was used to construct a linkage map [56]. Markers that were identical across all genotypes were identified and eliminated as duplicates. Furthermore, χ2-test for goodness-of-fit (*p* ≤ 0.0001) was used to identify markers with distorted segregation patterns [50,57,58]. Markers with significant deviation from the expected Mendelian segregation ratio (1:2:1) for F2:3 population were excluded from the analysis, resulting in a total of 1139 SNP markers used for the genetic map construction and QTL analysis.

#### *2.4. QTL Analysis and Candidate Gene Identification*