Development of a High-Density SNP-Based Linkage Map and Detection of QTL for β-Glucans, Protein Content, Grain Yield per Spike and Heading Time in Durum Wheat

High-density genetic linkage maps of crop species are particularly useful in detecting qualitative and quantitative trait loci for important agronomic traits and in improving the power of classical approaches to identify candidate genes. The aim of this study was to develop a high-density genetic linkage map in a durum wheat recombinant inbred lines population (RIL) derived from two elite wheat cultivars and to identify, characterize and correlate Quantitative Trait Loci (QTL) for β-glucan, protein content, grain yield per spike and heading time. A dense map constructed by genotyping the RIL population with the wheat 90K iSelect array included 5444 single nucleotide polymorphism (SNP) markers distributed in 36 linkage groups. Data for β-glucan and protein content, grain yield per spike and heading time were obtained from replicated trials conducted at two locations in southern Italy. A total of 19 QTL were detected in different chromosome regions. In particular, three QTL for β-glucan content were detected on chromosomes 2A and 2B (two loci); eight QTL controlling grain protein content were detected on chromosomes 1B, 2B, 3B (two loci), 4A, 5A, 7A and 7B; seven QTL for grain yield per spike were identified on chromosomes 1A, 2B, 3A (two loci), 3B (two loci) and 6B; and one marker-trait association was detected on chromosome 2A for heading time. The last was co-located with a β-glucan QTL, and the two QTL appeared to be negatively correlated. A genome scan for genomic regions controlling the traits and SNP annotated sequences identified five putative candidate genes involved in different biosynthesis pathways (β-glucosidase, GLU1a; APETALA2, TaAP2; gigantea 3, TaGI3; 14-3-3 protein, Ta14A; and photoperiod sensitivity, Ppd-A1). This study provides additional information on QTL for important agronomic traits that could be useful for marker-assisted breeding to obtain new genotypes with commercial and nutritional relevance.


Introduction
Durum wheat (Triticum turgidum L. subsp. Durum) is the tenth most widely grown crop on the planet with production of more than 36 million metric tons. It is commonly cultivated on both south and north shores of the Mediterranean Sea, with the largest proportion of the global production

Field Trait Analysis
Genetic variation in β-glucan content (BG), grain protein content (GPC), grain yield per spike (GYS) and heading time (HT) was assessed in a Duilio × Avonlea (D × A) RIL population grown at Policoro (Matera, Italy) and Valenzano (Bari, Italy) in 2014 ( Figure 1). All traits were determined in three field replicates in order to evaluate the environmental influence. Descriptive statistics (general mean for each parent and RIL, range of the RIL population, genetic variance, and broad-sense heritability) are reported in Table 1. Table 1. Descriptive statistics of BG, GYS, GPC and HT for the two parental lines (Avonlea and Duilio) and the RIL population D × A grown at Policoro (Matera) and Valenzano (Bari) in 2014. The genetic variance (s 2 G) and broad-sense heritability (h 2 ) of the traits are reported for the RIL population.  All traits were determined in three field replicates in order to evaluate the environmental influence. Descriptive statistics (general mean for each parent and RIL, range of the RIL population, genetic variance, and broad-sense heritability) are reported in Table 1. The phenotypic data for yield and adaptive components highlighted considerable differences between the parents and significant genetic variation in the RIL population in both environments. Avonlea was later heading and had lower GYS than Duilio. Qualitative trait data showed differences in BG, with cv. Avonlea having higher content, whereas GPC was similar for both parents. For all traits, the range of the RIL population was much larger than the range of the parental lines, suggesting transgressive segregation and that favourable alleles were present in both parents.
The means of the RIL population for BG were near Duilio, but with a big range of variability (0.30-0.63% and 0.22-0.68%) in both the environments with a typical pattern of a quantitative trait. The broad-sense heritability estimates (genotype mean basis) ranged from 0.80 to 0.82, indicating the stability of the trait and the phenotypic expression mainly due to genotypic effects.
Protein content in the RIL population ranged from 10.2 to 17.6% with mean values corresponding to the parental values. Broad-sense heritability (0.52 and 0.44) indicated the degree to which phenotypic expression was influenced by genetic versus environmental effects.
Similarly low estimates of heritability were observed for GYS. Moreover, GYS data showed a wide range of variability in the D × A RIL population (0.94-2.64 g and 1.58-3.47 g), typical of a quantitative trait and with mean values corresponding to the parental means. There was wide variation and transgressive segregation for HT with a range of values from 9 to 41. The heritability was 0.99-0.98, indicative of a high degree of stability for this trait. Contrary to the other traits analysed, heading date was bimodally distributed.
Correlation analyses (Table 2) between the various traits indicated significant and positive relationships for GPC and HT in Valenzano 2014 (p ≤ 0.01). A significant and negative correlation was detected between GPC and BG (p ≤ 0.01) in Valenzano 2014. Negative correlations between GPC and GYS were observed at both locations (p ≤ 0.001).

A Duilio × Avonlea Linkage Map
The D × A RIL population was genotyped with the wheat 90K SNP iSelect assay [7]. Of 81,587 assays, 5577 (6.8%) failed and 69,692 (85.4%) were monomorphic across the mapping population. The remaining 6831 (8.37%) were polymorphic; however, 513 had more than 10% missing data and were excluded from further analyses. Hence, 6318 markers were used for map construction.
Of 6318 loci, 5444 (6.7% from the original SNP dataset) assembled in 36 linkage groups when using a LOD score ≥ 6, and 874 remained unlinked or assembled in small linkage groups (LGs) without a chromosomal anchor locus. Linkage groups were assigned to the A and B genome chromosomes using the physical position of SNP reported in Colasuonno et al. [13] and the SNP anchor loci reported in the durum consensus map of Maccaferri et al. [15] (Figure 2, Table S1). Of 6318 loci, 5444 (6.7% from the original SNP dataset) assembled in 36 linkage groups when using a LOD score ≥ 6, and 874 remained unlinked or assembled in small linkage groups (LGs) without a chromosomal anchor locus. Linkage groups were assigned to the A and B genome chromosomes using the physical position of SNP reported in Colasuonno et al. [13] and the SNP anchor loci reported in the durum consensus map of Maccaferri et al. [15] (Figure 2, Table S1).    (Table 3). Linkage analysis showed a high number of co-segregating SNP (1596). The SNP markers were generally well distributed throughout the genome, although some chromosomes exhibited higher densities.
The overall SNP density was 2.7 markers/cM, with a maximum of 4.6 for chromosome 7B and a minimum of 1.6 for chromosome 3B. The recombinant frequency was 49.7%. The percentage of mapped SNPs showing segregation distortion at p ≥ 0.05 was 5.6%.

Detection of QTL for β-Glucan Contents, Protein Content, Grain Yield per Spike and Heading Time
Identification of QTL associated with BG, GPC, GYS and HT was carried out based on the genetic map obtained as previously described. Composite interval mapping (CIM) as proposed by Zeng [18] identified three QTL for β-glucans with LOD values ≥3.0 on chromosomes 2A and 2B, consistent in the mean of the two environments and in one location (Policoro), respectively ( Figure 2, Table 4). For additive QTL effects, positive and negative signs indicate the relative trait contributions from Avonlea and Duilio, respectively. The QTL on chromosome 2A had a negative effect due to the Duilio alleles, whereas both QTL on chromosome 2B showed positive effects by Avonlea alleles. Table 4. Marker-trait associations for BG, GPC, GYS, and HT determined for each location (Policoro and Valenzano) and for the mean values of the two environments. QTL, closest marker, chromosome locations, interval (cM), marker effects (additive effect with positive and negative value indicating contribution of parent Avonlea and Duilio, respectively), LOD value (logarithm-of-odds) ≥3.0, and R 2 (% of phenotypic variation explained by QTL) are reported for each SNP associated with each trait.  Eight QTL for GPC were located on chromosomes 1B, 2B, 3B (two loci), 4A, 5A, 7A and 7B. Positive alleles in Avonlea were on chromosomes 1B, 3B (both loci) and 7B, and those on 2B, 4A and 5A were from Duilio. The 7B QTL (QGpc.mgb-7B.1) was detected in both two environments and the mean.
Regression analysis allowed identification of a single region on chromosome 2A associated with HT that was consistent in both environments with a high LOD score of 25.5.

Candidate Genes Related to All the Detected QTL
Sequences of SNP markers in Wang et al. [7] were used to identify potential candidate genes associated with QTL for BG, GPC, GYS and HT detected in the D × A population. All the sequences of the SNP markers located in the QTL regions were analysed for annotated genes or proteins involved in each considered pathway in wheat. For BG, the marker IWB23783 located on chromosome 2B corresponded to β-glucosidase (GLU1a) ( Table 5). For GYS, candidate genes were detected on chromosomes 2B, 3A and 3B. Since the low number of SNPs in the 2B QTL region limited this study, we considered the same region in the consensus durum map [15]. This allowed us to detect the APETALA2 (TaAP2) gene corresponding to the wmc213-2B and wmc243-2B markers. Both SSRs mapped at 4.0 cM from the closest associated (IWA8152) marker. The gigantea 3 (TaGI3) gene showed sequence similarity with the IWB65703 marker (contig_3AS_200_3324175) at 0.6 cM from IWB48828 in the 3A QTL region. The 14-3-3 protein (Ta14A) gene, identified by the SNP marker IWB66165, mapped in the 3B QTL region at 7.0 cM from the closest SNP marker IWB64877 ( Table 5).
The heading time gene QHt.mgb-2A.1 appeared to be associated with the photoperiod sensitivity (Ppd-A) gene through SNP marker IWB54033 in the chromosome 2A QTL region.

Discussion
The recent emphasis on grain quality over grain quantity has shifted the major objective of wheat breeding programs to increasing nutrient properties while maintaining or increasing grain yield in lines to be released for commercial production. For a long time the major quality objective in durum breeding has been increased GPC [19] and yellow pigment content [9], but recently, more attention has been given to fibre [11], for which the main objective has been improved BG [10]. Different components contribute to grain yield, and among them the fundamental one is grain yield per spike, which is generally negatively correlated to GPC [2]. Another trait correlated with yield is HT, which in cereal crops is determined by vernalisation requirement, photoperiod sensitivity, and narrow-sense earliness or earliness per se [20]. Thanks to the recently released consensus linkage maps for bread wheat [8] and durum [15], development and saturation of genetic maps for detection of QTL underlying agronomic traits has rapidly increased in genetics studies.
In the present work, we evaluated BG, GPC, GYS and HT in an Italian cv. Duilio × Canadian cv. Avonlea RIL population grown in two field trials and developed a high-density linkage map. The total map length was 1962.5 cM, and included 5444 SNP markers, of which 374 SNPs were used to anchor the linkage groups (LGs) to chromosomes. A good coverage and representation of all 14 durum chromosomes was obtained, but the presence of gaps (over 20 cM) led to 36 LGs. Despite the high density of polymorphic markers detected using the 90K SNP assay, linkage groups ranging from 20 to 33 cM have been reported in previous SNP maps [13,[21][22][23]. This fragmentation observed for the D × A map could be dependent upon the strict conditions used in determining the marker order (LOD score > 6 and SNP markers with less than 10% missing data). Despite the numerous LGs, our map represented a useful basis for studies on marker-assisted selection and positional cloning. The D × A map had a higher length (1962.5 cM) than previous maps [2,3,13], but a lower marker density (2.7 cM). The average number of mapped markers per chromosome was 389, with a range from 230 SNPs on chromosome 4B to 536 on chromosome 7B. More markers were mapped on the B genome as observed in other linkage studies [2,3,13]. Chromosome 1A was represented by only one LG [15].
A total of 19 QTL for the examined traits were detected, including three for BG, eight for GPC, seven for GYS, and one for HT.
BG QTL were mapped on chromosomes 2A and 2B, as previously reported by Marcotuli et al. [10], but in different chromosome positions. QGbg.mgb-2A.1 was co-located with QHt.mgb-2A.1 controlling HT and with a strong and negative effect on the trait (−9.96) due to Duilio contributing the negative effect. Previous work on yield components reported a major QTL for HT on chromosomes 2BL with negative alleles on the trait [24,25]. We found no significant correlation between BG and HT. Similar results were reported in oat [26,27] and barley [28,29]. Due to the paucity of studies on BG in wheat at any ploidy level, the present data provide some knowledge and understanding of the trait in durum wheat.
To identify coding sequences in regions associated with QTL, BLAST analyses on the NCBI database [48] were carried out. No gene co-located with the BG QTL QGbg.mgb-2A.1 and QGbg.mgb-2B.1, whereas the mapped SNP IWB23783 associated with QGbg.mgb-2B.2 corresponded to β-glucosidase 1a (GLU1a), a member of the glycoside hydrolase family 1 (GH1). This enzyme catalyses hydrolysis of the glycosidic linkage in glycosides, leading to formation of a hemiacetal or hemiketal sugar and a corresponding free aglycon [49]. We also detected the photoperiod sensitivity (Ppd-A1) gene in the QTL QHt.mgb-2A.1 region. The contig sequence of IWB54033 located in the 2A QTL corresponded to the Ppd-1A sequence, confirming previous results on the relationship between photoperiod response and HT [50,51]. No marker located in gene sequences for GPC was identified, but previous studies reported genes involved in the nitrogen pathway located in the same position of QTL reported in the present work [2,[52][53][54][55].
Bioinformatics analysis of markers located in QTL region QGys.mgb-2B.1 for GYS detected no candidate gene associated with the trait. Because previous studies reported a gene (APETALA2) associated with a QTL for grain yield in different species [56][57][58][59], we decided to use all the markers (SNPs and Simple Sequence Repeat, SSR) from the durum wheat consensus map [15] co-locating with all the SNP mapped in the 2B chromosome interval 19.9-37.9 cM of the D × A map. Due to the higher number of considered markers, we identified the APETALA2 (TaAP2) gene associated with the SNP marker IWA8152.

Genetic Materials and Field Experiments
A set of 134 F 2:7 RILs were developed from a cross between Italian cultivar Duilio and Canadian cultivar Avonlea. After the last generation of selfing, each line was harvested in bulk to provide seeds for replicated trials and DNA extraction. The parental lines and RILs were evaluated for β-glucan (BG), grain protein content (GPC), grain yield per spike (GYS), and heading time (HT) in replicated trials at Policoro (metropolitan city of Matera) and Valenzano (metropolitan city of Bari) in southern Italy in 2014.
A randomised complete block design with three replications was used in field experiments. Plots consisted of 1 m rows spaced 30 cm apart; each plot had 80 germinated seedlings. During the growing season, 10 g of nitrogen per m 2 were applied at flowering stage and standard cultivation practices were adopted.
Plots were hand-harvested at maturity and GYS was determined by dividing grain yield by the number of spikes (about 70-80 spikes) per row.
β-glucan (percentage w/w) content in the whole grain was assayed using a Mixed-Linkage β-Glucan Assay Kit (Megazyme International Ireland Ltd, Wicklow, Ireland) based on the method by McCleary and Codd [62] and included the industrial standard for barley (4.1% β-glucan). GPC, expressed as a percentage of protein on a dry weight basis, was determined on a 2 g sample of whole-meal flour using an InfraAlyzer spectrophotometer (Technicon Industrial Systems, Tarrytown, NY, USA) (that is, using near-infrared reflectance spectroscopy).
HT was recorded as the number of days from 1 April 2014, to 50% ear-emergence, corresponding to stage 55 on the Zadoks et al., scale [63].

SNP Genotyping
Fifty ng/µL of genomic DNA from each RIL and parental line (Duilio and Avonlea) was analysed with the wheat 90K iSelect array [8]. Genotyping was performed at TraitGenetics GmbH (http://www.traitgenetics.de [64]) following the manufacturer's recommendations as described in Akhunov et al. [65]. The genotyping assays were carried out using the Illumina iScan reader and performed using GenomeStudio software version 2011.1 (Illumina, San Diego, CA, USA).

Segregation Analysis and Map Construction
Chi-squared tests were used to determine the goodness-of-fit at p > 0.001 of segregation ratios to expected 1:1 ratios for each SNP. All markers with more than 10% missing data, or segregating presence/absence in the mapping population were excluded from further analysis. Linkage analysis between markers and determination of the linear order of loci was performed by JoinMap 4.0 (Kyazma, Wageningen, Netherlands) [66] using the regression mapping algorithm. Grouping was performed using the independence LOD parameter, with groups presenting a LOD 6. The Haldane mapping function was used to calculate map distances [67]. Physically mapped SNP markers [13] and SNP data from the durum consensus maps [15] were used as anchor loci and for assigning linkage groups to specific chromosomes. Linkage groups were named according to the wheat chromosome nomenclature followed by a number linkage groups within chromosomes.

Statistical Analysis, and QTL and Genes Detections
Analysis of variance was carried out for each trait using GenStat (version 18, VSN International Ltd, Hemel Hempstead, UK,). QTL detection was performed in QGene 8.3.16 using composite interval mapping [18]. A marker-trait association was considered significant when one or more markers was associated with a trait at a threshold of −log 10 (p) ≥ 3.0, determined by modified Bonferroni correction. For additive effects, positive and negative signs for estimates indicated the contributions of Duilio and Avonlea, respectively, toward higher trait values. Regression analysis implemented in QGene was used for HT. Graphical representations of linkage groups and QTL were carried out using MapChart 2.2 software (SolarWinds, Austin, TX, USA) [68].
Putative candidate genes associated with QTL for BG, GPC, GYS and HT were identified using SNP sequences from Wang et al. [8] blasted against Triticum sequences annotated in NCBI [69] and contig sequences from the Unité de Recherche Génomique Info (URGI) website [70].

Conclusions
In the present study, we report a high-density SNP-based genetic map for the Duilio × Avonlea durum wheat RIL population. It comprised 5444 SNP loci. This high-resolution map and the data on β-glucan content, protein content, grain yield per spike, and heading time evaluated in two environments provided increased knowledge of BG in durum wheat and its interaction with other quantitative traits. The three QTL for BG were influenced by environmental differences. The QTL on chromosome 2B was coincident with the β-glucosidase 1a gene involved in the BG pathway. The QTL for BG on chromosome 2A was coincident with HT providing additional information on trait correlations. The marker-trait associations for GPC and GYS confirmed previous data reported in the literature on the polygenic nature of the traits and correspondent QTL in both common and durum wheat. The overall results provide a basis for marker-assisted selection for BG in breeding programs.