Identification of Quantitative Trait Loci Relating to Flowering Time, Flag Leaf and Awn Characteristics in a Novel Triticum dicoccum Mapping Population

Tetraploid landraces of wheat harbour genetic diversity that could be introgressed into modern bread wheat with the aid of marker-assisted selection to address the genetic diversity bottleneck in the breeding genepool. A novel bi-parental Triticum turgidum ssp. dicoccum Schrank mapping population was created from a cross between two landrace accessions differing for multiple physiological traits. The population was phenotyped for traits hypothesised to be proxies for characteristics associated with improved photosynthesis or drought tolerance, including flowering time, awn length, flag leaf length and width, and stomatal and trichome density. The mapping individuals and parents were genotyped with the 35K Wheat Breeders’ single nucleotide polymorphism (SNP) array. A genetic linkage map was constructed from 104 F4 individuals, consisting of 2066 SNPs with a total length of 3295 cM and an average spacing of 1.6 cM. Using the population, 10 quantitative trait loci (QTLs) for five traits were identified in two years of trials. Three consistent QTLs were identified over both trials for awn length, flowering time and flag leaf width, on chromosomes 4A, 7B and 5B, respectively. The awn length and flowering time QTLs correspond with the major loci Hd and Vrn-B3, respectively. The identified marker-trait associations could be developed for marker-assisted selection, to aid the introgression of diversity from a tetraploid source into modern wheat for potential physiological trait improvement.


Introduction
Wild relatives and ancestral landraces of wheat offer a plethora of genetic diversity that could be utilised in breeding for improved yield and environmental stability in modern wheat. Domestication from progenitor species, coupled with modern breeding techniques, has narrowed the genepool of recent varieties [1]. Domestication of wheat took place around 12,000 BP and included diploid wheat as well as domesticated forms of the tetraploid T. dicoccoides, predominately Triticum turgidum subsp. dicoccum (T. dicoccum), commonly known as emmer wheat [2]. Modern wheat grown today mainly belongs to two species: the hexaploid Triticum aestivum (2n, 6x, 42 chromosomes; for awn length (A L ) and median for trichome density (T D ). Only for stomatal density (S D ) did the variation between parents differ significantly across years. In 2017, dic12b had the significantly higher mean of 74.2 mm −2 . However, Tios had the higher mean of 72.4 mm −2 in 2019, although the difference was not significant (p = 0.09). Across both trials, dic12b had a higher mean FL L than Tios, although the difference was only significant in 2019. For A L , S D and T D , generalised heritability (H 2 ) was higher in the 2017 trial compared to 2019 (Table 1). In 2017, H 2 could not be calculated for FL L due to the negative variance component associated with genotype in the model. Across both years, A L was the most heritable trait. An example of A L diversity within the population and between parents is shown in Supplementary Figure S1. In 2017, FL W had the lowest calculated H 2 (0.36), whilst in 2019, the lowest heritability was observed for S D (0.36). Table 1. Trait means for each parent (Tios and dic12b) and the tetraploid mapping population (PS1) individuals across the 2017 field trial and the 2019 pot trial. Generalised heritability (H 2 ) is shown for each trait. Mean values in the table were calculated using raw data. Within each year, the asterisks by each trait name indicates significant differences at particular p-values between the mapping population parents (Tios and dic12b). The full test results are shown in Supplementary Tables S1 and S2.  Pearson's correlation coefficient values between traits across both trials are shown in Figure 1. A significant positive correlation (p < 0.01) was observed for every trait when measured across the different years. A L and F T showed the strongest correlations between years (r = 0.86 and 0.64, respectively), while FL L and FL W showed the weakest correlations between years (r = 0.41 and 0.43, respectively). FL L from 2019 was significantly negatively correlated with F T from the 2017 and 2019 trials, while FL L from 2017 was only negatively correlated with F T in 2017. In 2019, FL L also showed a significant negative correlation with S D , while a positive correlation was observed with FL W . Additionally, a negative correlation was found in 2019 between S D and T D . This relationship was not observed in 2017, where S D and T D were significantly positively correlated. Finally, a positive relationship was also observed between FL L and A L in the 2017 trial.

Linkage Map
Using 104 genotyped PS1 individuals, a new genetic linkage map was formed consisting of 2066 markers along 14 chromosomes (Figure 2a). Within each chromosome, the new genetic positions (cM) of the markers were correlated against the physical positions (Mb) for the same SNPs. Visual inspection of these plots indicated both appropriate ordering and orientation of each chromosome within the new genetic map. Examples of these plots are shown in Figure 2b for chromosomes 6A and 3B, which had the lowest and highest frequency of markers per chromosome in the linkage map (82 and 214 markers, respectively). Two heat maps were plotted as a further diagnostic of the new map: the pairwise recombination fraction and the logarithm-of-odds (LOD) linkage between ordered markers, shown in the top and bottom triangles of Figure 2c. Low recombination fractions and high LOD scores of closely neighbouring ordered SNPs indicated that the map was formed appropriately, this is shown through the high heat observed along the central diagonal of both heat maps in Figure 2c. . Figure 1. Pearson's correlation matrix for the 6 traits measured across both trials. Where a square is blank, the Pearson's product moment correlation test yielded a non-significant p-value (p > 0.01). AL = awn length; FT = flowering time; FLL = flag leaf length; FLW = flag leaf width; SD = flag leaf stomatal density; TD = flag leaf trichome density. The correlation matrix was formed in R using the package corrplot [31]. Plot font was edited using the R package extrafont [32].

Linkage Map
Using 104 genotyped PS1 individuals, a new genetic linkage map was formed consisting of 2066 markers along 14 chromosomes (Figure 2a). Within each chromosome, the new genetic positions (cM) of the markers were correlated against the physical positions (Mb) for the same SNPs. Visual inspection of these plots indicated both appropriate ordering and orientation of each chromosome within the new genetic map. Examples of these plots are shown in Figure 2b for chromosomes 6A and 3B, which had the lowest and highest frequency of markers per chromosome in the linkage map (82 and 214 markers, respectively). Two heat maps were plotted as a further diagnostic of the new map: the pairwise recombination fraction and the logarithm-of-odds (LOD) linkage between ordered markers, shown in the top and bottom triangles of Figure 2c. Low recombination fractions and high LOD scores of closely neighbouring ordered SNPs indicated that the map was formed appropriately, this is shown through the high heat observed along the central diagonal of both heat maps in Figure  2c.
The total length of the new genetic map spanned 3295 cM, with an average spacing (between markers) of 1.6 cM. The largest spacing of 26.6 cM was found on chromosome 6A (Figure 2a). Chromosome length varied from 147.8 cM (6A) to 303.8 cM (7A). The lowest density of markers was found on chromosome 4A, with an average spacing (between markers) of 2.3 cM, and the highest density on chromosome 1B, with an average spacing of 1.1 cM. Across the new map, the B genome had a higher coverage of markers than the A genome, with 1185 associated markers across an average spacing of 1.4 cM compared to 881 markers with an average spacing of 1.9 cM. The complete genetic map formed using the PS1 population is shown in Supplementary Table S3.

Figure 1.
Pearson's correlation matrix for the 6 traits measured across both trials. Where a square is blank, the Pearson's product moment correlation test yielded a non-significant p-value (p > 0.01). A L = awn length; F T = flowering time; FL L = flag leaf length; FL W = flag leaf width; S D = flag leaf stomatal density; T D = flag leaf trichome density. The correlation matrix was formed in R using the package corrplot [31]. Plot font was edited using the R package extrafont [32].
The total length of the new genetic map spanned 3295 cM, with an average spacing (between markers) of 1.6 cM. The largest spacing of 26.6 cM was found on chromosome 6A (Figure 2a). Chromosome length varied from 147.8 cM (6A) to 303.8 cM (7A). The lowest density of markers was found on chromosome 4A, with an average spacing (between markers) of 2.3 cM, and the highest density on chromosome 1B, with an average spacing of 1.1 cM. Across the new map, the B genome had a higher coverage of markers than the A genome, with 1185 associated markers across an average spacing of 1.4 cM compared to 881 markers with an average spacing of 1.9 cM. The complete genetic map formed using the PS1 population is shown in Supplementary Table S3.

QTL Mapping
Results from the QTL mapping for both trials are shown in Table 2. In the 2017 trial, six candidate QTLs were identified for four traits across five chromosomes (2B, 4A, 5B, 6A and 7B). The percentage phenotypic variation explained by each QTL ranged from 12.3% (Q3_A L _17) to 35.1% (Q5_A L _17), and LOD scores ranged from 4.4 (Q8_FL L _17) to 12.2 (Q5_A L _17). The second highest LOD score was observed for Q9_F T _17 (LOD = 7.1). The size of the 1-LOD support intervals, expanded to the nearest flanking marker, ranged from 7.2 cM for Q7_FL W _17 to 124.5 cM for Q3_A L _17. Of the QTLs found in the initial single genome scan for each trait, only Q7_FL W _17 fell below the 0.05 alpha LOD threshold, but still passed the 0.1 threshold with an initial LOD score of 4.2. The additive effect for each QTL, shown in Table 2, indicated the effect of one dic12b allele (0.5 of BB) to the trait data. There was a positive additive effect of a dic12b allele at the QTLs: Q3_A L _17, Q5_A L _17 and Q7_FL W _17, whilst there was a negative additive effect of a dic12b allele at the QTLs: Q2_FL W _17, Q8_FL L _17 and Q9_F T _17. Using the residuals of F T regressed on the FL L data from 2019 produced no significant QTLs and, using the same approach with the 2017 data, there was little difference from the existing model found for Q9_FT_17. Using the residuals of the FL L data regressed on F T from 2019 produced no QTLs. However, using the residuals of the FL L data regressed on F T from 2017 produced a QTL on 6A (Q8_FL L _17) that was not found using the raw trait data. The QTL had a LOD score marginally above the 5% alpha of 4.34 (LOD = 4.4). The parent dic12b had greater FL L across both years, although a dic12b allele at Q8_FL L _17 caused a negative effect on the trait (Table 2). However, the variation in FL L between parents was reduced in the 2017 trial compared to 2019.  [35] for two chromosomes: 3B and 6A, which showed the highest and lowest marker frequency, respectively. (c) A heat map showing the assembled PS1 linkage map. The logarithm-ofodds (LOD) linkage between ordered markers is plotted in the bottom corner and pairwise recombination fractions between ordered markers is shown in the upper triangle, the heat map was created using the 'heatMap' function in ASMap [33]. Plot fonts were edited using the R package extrafont [32]. The logarithm-of-odds (LOD) linkage between ordered markers is plotted in the bottom corner and pairwise recombination fractions between ordered markers is shown in the upper triangle, the heat map was created using the 'heatMap' function in ASMap [33]. Plot fonts were edited using the R package extrafont [32]. Table 2. Candidate quantitative trait loci (QTLs) identified using the tetraploid mapping population (PS1) for five traits, using phenotype data from two separate trials. Logarithm-of-odds (LOD) scores, percentage phenotypic variation explained by QTLs (% var) and additive effects were extracted from models fitted through the 'fitqtl' function in R/QTL. The interval shown for each QTL is the 1-LOD support interval expanded to the flanking markers. LOD thresholds are shown at the significance levels of 5% Alpha and 10% Alpha. These were determined through running 5000 permutations within the R/QTL 'scanone' function implemented through Haley-Knott regression [36]. * Second and third elements of each QTL name signified the related trait and trial, respectively. ** Additive effect of 0.5 BB genotype (dic12b parent). + Data were expressed to the log10. ++ Residuals were used as trait data: the F T 2017 data were corrected for germination heterogeneity and the FL L 2017 data were corrected for F T .

QTL
In the 2019 trial, four candidate QTLs were identified for four different traits, located on four chromosomes (Table 2). In 2019, Q4_A L _19 yielded the largest LOD score and phenotypic variation explained by a QTL (11.6 and 39%, respectively). The distribution of the QTLs mapped across the genome in both years is shown in Figure 3. Several QTLs collocated: Q4_A L _19 was mapped to 11 cM on 4A whilst Q5_A L _17 mapped to 15.1 cM on 4A in the 2017 trial. Q10_F T _19 was mapped to 20 cM on 7B, which co-locates with Q9_F T _17 in the 2017 trial which was mapped to 13 cM on 7B. Further evidence that consistent QTLs were found across multiple years for A L and F T was shown through the dic12b allele contributing a positive additive effect at both Q4_A L _19 and Q5_A L _17, while the dic12b allele at both Q10_F T _19 and Q9_F T _17 contributed to a decrease in F T . Additionally, Q6_FL W _19 mapped to 206.7 cM on 5B, was likely to be the same QTL as Q7_FL W _17 mapped to 226.5 cM in 2017, where the dic12b allele contributed a positive effect on the trait for both QTLs. A single borderline QTL was found for T D in 2019 (Q1_T D _19) with a LOD score of 4.4, which was only slightly higher than the LOD 5% alpha threshold of 4.34.
Plants 2020, 9, x FOR PEER REVIEW 7 of 17 dic12b allele at both Q10_FT_19 and Q9_FT_17 contributed to a decrease in FT. Additionally, Q6_FLW_19 mapped to 206.7 cM on 5B, was likely to be the same QTL as Q7_FLW_17 mapped to 226.5 cM in 2017, where the dic12b allele contributed a positive effect on the trait for both QTLs. A single borderline QTL was found for TD in 2019 (Q1_TD_19) with a LOD score of 4.4, which was only slightly higher than the LOD 5% alpha threshold of 4.34.  The physical positions of the peak SNPs and 1-LOD support interval flanking markers for each QTL are listed in Table 3. As with genetic positions (Table 2), the largest physical 1-LOD support interval was found for Q3_A L _17 (653 Mb), whereas the smallest physical interval (17.29 Mb) was found for Q9_F T _17 ( Table 3). The peak physical marker positions of the QTLs for A L mapped to 4A across different years (Q4_A L _19 and Q5_A L _17) were relatively distant. However, the physical 1-LOD support intervals for each QTL were similar: 40 to 103 Mb for Q4_A L _19 and 47 to 110 Mb for Q5_A L _17. The peak marker was identical for Q9_F T _17 and Q10_F T _19. However, the peak markers for the QTLs mapped across different years for FL W on 5B were 52 Mb apart. Table 3. Genetic and physical positions of peak single nucleotide polymorphisms (SNPs) closest to quantitative trait loci (QTLs) identified using the tetraploid mapping population (PS1). Physical positions are also shown for the 1-LOD support interval markers. An NA indicated where an interval marker aligned to a different chromosome than the one mapped in the genetic map. Physical positions were taken from Basic Local Alignment Search Tool (BLAST+) results against RefSeq v1.0 wheat genome assembly [35].

Discussion
Typically, a polygenic trait will follow a relatively normal distribution, which was the case for all traits except trichome density (T D ). In bread wheat, awn development is considered to be determined by a few major loci [38], while the flowering time pathway is considered to be more complex [39]. However, transgressive segregation appeared greater in awn length (A L ) than flowering time (F T ) across both trials (Supplementary Figure S2). Transgressive segregation was still observed for every trait, indicating the presence of complementary alleles inherited in the progeny from both parents [40]. This combination of complementary alleles, leading to progeny with extreme phenotypes, would be useful for trait-specific ideotype selection, and hence genetic improvement of these traits in breeding programmes. Furthermore, high heritability was found consistently for three of the six traits measured (A L , F T and T D ) and moderate heritability was observed for the remaining traits (stomatal density (S D ), flag leaf length (FL L ) and width (FL W )) in at least one trial. The high to moderate heritability observed for the traits indicated high genetic variability within the population and suggested that there were likely to be underlying QTLs linked to phenotypic variation (particularly for A L , F T and T D ).
For A L and F T , there were major loci segregating in the population, which meant that these QTLs were easily detectable. For both S D and FL L, there was no consistent significant variation observed between the parent lines across trials. Furthermore, low H 2 was observed for S D in 2019 and no H 2 estimate could be obtained in 2017 for FL L , due to a negative variance term associated with genotype. These observations signified that no major loci were segregating for these traits and that the variation observed may be a result of the environment. However, across the two trials, each of these traits were significantly correlated, suggesting the presence of some genotypic effects controlled by minor QTLs. This was supported by the 2017 trial, where there were moderate H 2 observed for S D (Table 1) and a minor QTL found for FL L (Table 3). It was probable that the relatively small population size used reduced the power of detection of further minor QTLs. The distribution of the T D data before the log transformation suggested that there may have been a major loci controlling the trait, further supported by the consistent significant variation between parents across trials and the high H 2 observed. However, only a single minor QTL was found for T D in 2019 that contributed a small effect (17%) to the observed phenotypic variation, potentially due to poor marker coverage surrounding the segregating loci. There have been limited investigations into QTLs controlling T D in wheat. The negative correlation observed between T D and S D in 2019 indicated that the T D variation may have been influenced by physiological trade-offs with other phenotypes.
One of the difficulties in ideotype design and trait-based selection are the potential trade-offs or compensations between characteristics [18,19]. This was apparent in this study from a number of negative correlations observed in the trials (Figure 1), such as between FL L and F T , which was observed over multiple years. No QTLs were found for S D in either trial. In 2019, the Tios parent had a noticeably shorter flag leaf than in 2017, and as FL L negatively correlated with S D , this may have caused an increase in S D on the flag leaves of the Tios line. T D correlated negatively with S D in 2019, which could suggest that the partitioning of more assimilates towards a higher stomatal density per unit leaf area may cause a reduction in trichome number. This negative relationship could be beneficial for breeding in a water-rich environment, where stomatal density would be favoured over trichome density for increased CO 2 capture, regardless of water loss. However, a significant positive correlation was observed between the traits in 2017 (Figure 1), suggesting that environmental factors may have had the greatest influence on the relationship. There were trait relationships reflecting assimilate investment into both the flag leaf and the ear: FL L and A L were positively correlated in the 2017 trial. Investment into ear area has been previously linked to reductions in flag leaf area [41], so the combination of increased organ size found in the present study may be promising for breeding high yielding lines [42]. Marker-trait associations, such as the ones found in the present study, would be useful tools in the selection for maintaining favourable trait combinations.
Three consistent QTLs were found over multiple years linked to three traits (A L , F T and FL W ). The most significant of these was linked to A L and mapped to 4A, where the longer awned parent (dic12b) allele increased A L ( Table 2). Although not yet cloned, there is a known gene on the short arm of 4A that causes awn shortening and bending in common wheat, called the Hd locus [43]. In a hexaploid wheat recombinant inbred line population, Yoshioka et al. [38] identified three major QTLs related to awn length which corresponded to the three known inhibitor genes: Hd, B1 and B2 (located on 4A, 5A and 6B, respectively). Yoshioka et al. [38] found that the QTL on the short arm of 4A corresponded to the Hd locus located at 22 cM and explained 23.1% of the phenotypic variation. In another hexaploid wheat population, Sourdille et al. [44] mapped two QTLs relating to awn length, one of which segregated with the Hd gene. The genetic control of awn length and presence is less studied in T. dicoccum. However, in a T. dicoccum and T. durum mapping population, Buerstmayr et al. [45] found a QTL on the short arm of chromosome 4A that was linked to A L and they concluded that it may correspond to the Hd locus. Although the confidence interval identified in the present study was relatively large (Tables 2 and 3), it is probable that the 4A QTL was segregating with the major Hd locus. The second QTL linked to A L was found on chromosome 2B in the 2017 trial, the same QTL was also apparent in 2019 but fell just under the significance thresholds used and was thus not reported. Other studies have also found QTLs controlling A L aside from the three known inhibitor gene loci [46,47]. However, there is no previous evidence in the literature supporting the 2B QTL found in the present study.
The F T QTL mapped to 7B had the smallest physical 1-LOD support interval of the QTLs mapped across multiple years (Table 3). Flowering time control in tetraploid wheat is less studied than in hexaploid wheat [48]. In hexaploid wheat, the genetic control of F T is complex and controlled by three signalling pathways: vernalisation (Vrn genes), photoperiod response (Pdp) and earliness per se (Eps) [39,49]. Würschum et al. [48] studied 328 European durum genotypes and found six QTLs linked to heading time, with their hits on 2B and 7B corresponding to the known genes Ppd-B1 and Vrn-B3, respectively. The physical position of the 7B QTL in Würschum et al. [48] fell within the physical 1-LOD support interval of the QTL mapped on 7B in the present study. Variation in Vrn-B3 has been found to be prolific in durum wheat [50] and in T. dicoccum [51]. Vrn-B3 is considered to be linked exactly to the gene in Arabidopsis Flowering Locus T (FT) and is known to have some control on flowering time [52]. Therefore, although both the mapping parents (Tios and dic12b) had no vernalisation requirement, it is probable that this major locus was linked to flowering time in the PS1 population.
A QTL was only identified for FL L in 2017 once the effect of F T had been removed. As no additional QTLs were identified for F T once the FL L effect had been removed, it was probable that an earlier flowering time influenced flag leaf growth. F T is known to be an important factor in determining many physiological traits [53]. A stable QTL was found for FL W towards the end of 5B (Table 2). Other studies have identified a considerable number of environmentally stable QTLs for FL W and FL L in hexaploid wheat [54,55]. However, there has been limited investigation into QTLs linked to flag leaf area in a tetraploid wheat background. Yang et al. [56] did identify a QTL linked to FL W in a hexaploid population on the long arm of 5B which was stable under multiple environments, although they mapped 15 QTLs for FL W in total. This evidence indicates that FL W in wheat is probably controlled by many minor QTLs, including the 5B QTL identified in the present study.
The first step of MAS application in breeding programs is the identification of candidate QTLs for marker-trait associations, which has been achieved in the present study for five hypothesised proxies that relate to photoassimilate supply or drought tolerance. Compared to modern wheat, there has been relatively limited previous investigation into the genetic control of this combination of traits in tetraploid landraces. The genotyped T. dicoccum mapping population produced in this study is a useful tool for the research and breeding community and could be employed further for the identification of marker-trait associations in a diverse background. Improved resolution of QTL locations in this study would be needed before these markers could facilitate the introgression of diverse physiological characteristics into modern wheat from a tetraploid source. Building on the present study, this goal could be achieved by fine mapping approaches or by increasing the number of individuals in the mapping population.

Plant Material
A tetraploid mapping population (PS1) was created by crossing two T. dicoccum landraces, shown in Supplementary Figure S1. The landraces were taken from a collection described in Leigh et al. [57], and the maternal parent 'Tios' and the paternal parent 'dic12b' originated from Spain [58] and Iran, respectively. Parental samples were taken from single founder plants to reduce inherent genetic heterogeneity in these landraces. F2 seed of 200 individuals originating from six F1 plants was sown into 96-well trays and grown in glasshouses at Park Farm, Cambridge. Individuals were taken through a single seed decent program to the F5 generation from October 2014 to March 2017. At each generation, flowering ears were covered with a clear cellophane bag to prevent outcrossing and promote self-fertilisation. By the F5 generation, the population size had reduced to 101 due to germination failures and sterility, which may have been related to the poor adaption of the landraces to the glasshouse environment. As a result, F4 seed was included in a field trial with the F5 seed to increase the population size. Of the 120 individuals that reached maturity in the trial, 87 and 33 were from F5 and F4 seed, respectively. Seed was hand threshed from ears collected from this trial and 108 individuals were advanced to a further generation (F5 or F6) in an outdoor 5 L pot trial sown on 29 April 2019.

Trials and Phenotyping
The PS1 population was phenotyped in two trials conducted in outside areas. Firstly, a mix of F4 and F5 plants were grown in a field trial at the New Ornamentals Field, Park Farm, NIAB, Cambridge. Seeds were hand-sown into two replicated 25.5 m 2 plots (8.5 × 3 m) on 24 March 2017. Each plot had nine 8.5 m long rows spaced every 30 cm. The 1st and 9th row and the 25 cm at the start and end of each row was sown with seed of hexaploid bread wheat variety Paragon to create a buffer and to reduce edge effects. The 2nd, 5th and 8th rows were not sown with any seed, in order to aid accessibility for phenotyping. The 3rd, 4th, 6th and 7th rows were organised into a randomised incomplete block design with 5 blocks per plot. Within each block, there were 32 mini-plots (30 mapping individuals and the 2 parents). Each mini-plot was a 20 cm row segment, consisting of a label and 3 seeds sown 5 cm apart. The trial was fertilised with one application of ammonium nitrate (34.5% N) at a rate 434 kg ha −1 , in order to achieve a universal N application of 150 kg/N, and treated with insecticide and fungicide in two applications. The trial was manually irrigated. Every plant in the trial was individually staked to avoid lodging. Non-invasive measurements were made on the flag leaves of the primary tiller from the first individual in each mini-plot to reach growth stage (GS) 61 [59]. Leaf impressions were taken from the second plant to reach GS61. If no second plant germinated across the trial for a particular line, the impression was taken on the flag leaf of a secondary main tiller of the first plant. Only one of the two replicated plots was used for the measurements; however, if an individual failed to grow in the 1st plot it was measured from the 2nd plot. All parent mini-plots were measured across the trial in order to test for potential differences between plots.
Seed was harvested from the 2017 field trial and sown in a subsequent 2019 trial. Six seeds of the same individual were sown in one 5 litre pot on 4 April 2019 in an unheated glasshouse and thinned on germination to leave three viable seedlings. Pots were transferred to an outside trial area at Park Farm, NIAB, Cambridge, on 24 April 2019, where they were laid out in an incomplete block design with two replicate pots per individual, 20 replicates of each parent and a block size of eight pots. The design consisted of 15 columns, with the outer columns containing pots filled with Paragon plants to act as a buffer. The 2nd, 5th, 8th, 11th and 14th columns were left as 0.5 m walkways to aid accessibly for phenotyping. All plants were grown in an 80% peat-reduced and 20% bark compost mix, with a slow release fertiliser incorporated into the mix. A single fungicide application was made. Within each pot, measurements were taken from primary tillers of the plant furthest from the walkway at GS61. All pots were connected to an automated drip irrigation system ensuring even water availability across the trial.
For both trials, measurements of flag leaf length (FL L ), width (FL W ), awn length (A L ) and the number of days from sowing to flowering time (GS61, F T ) were recorded. An abaxial leaf impression from halfway up the flag leaf was made onto a microscope slide using Loctite Power Flex Ehyl-2-Cyanocrylate (Loctite, Düsseldorf, Germany). A uniform mesophyll region from each slide was imaged using a Leica DM2500 Optical Microscope (Leica Biosystems, Nussloch, Germany) at 5x magnification. Stomata and trichomes were counted from the 5x images in three different 0.3 mm 2 overlaid grid squares and the mean values were taken, enabling stomatal and trichome density to be estimated per mm 2 (S D and T D , respectively). Only fully developed trichomes were included in the counts, where the trichome appendage was clearly visible.

Trial Analysis
Many of the seeds across the mini-plots in the 2017 field trial failed to germinate. However, for all traits except F T , there was no evidence that the germination heterogeneity influenced the trait (data not shown). For F T , there was strong evidence that the germination heterogeneity influenced the trait. F T trait data were therefore regressed on the germination scores and residuals were extracted from the model. The residuals were used as the adjusted trait data for QTL mapping. As the majority of the data in the 2017 field trial came from single replicates, no adjustment for spatial trends were made and raw data were used in the QTL mapping for the five remaining traits.
The increased use of replicates in 2019 enabled all traits to be adjusted for spatial effects, as implemented via the R package SpATS [60]. After adjusting spatial trends through two-dimensional P-spline modelling, best linear unbiased estimations (BLUE) per line were extracted and used as phenotypic data in the QTL mapping. Frequency histograms were plotted in R with the package ggplot2 [61] and showed that T D data from both years were not normally distributed, as the most frequent observed T D was close to 0. To avoid violating underlying assumptions made through statistical testing, T D data from both years were expressed to the log10. A negative correlation was found between FL L and F T in both trials. To assess marker-trait associations independent of this relationship, each trait was regressed on the other within each year. The residuals were extracted from each model and tested as phenotypic data in the QTL mapping, and the results were compared using the uncorrected trait values.
For the 2019 trial, generalised heritability (H 2 ) was calculated using the 'getHeritability' function within SpATS, which implements a calculation suggested by Oakey et al. [62]. For the 2017 trial, H 2 was also calculated on a line mean basis, implemented by the VHERITABILITY procedure in Genstat [63], which employs a calculation of H 2 from Cullis et al. [64].

Genotyping
DNA was extracted from seven-day-old seedlings of 115 F4 PS1 individuals and parental lines following a modified Tanksley 96-well-microprep extraction method [65]. DNA samples were sent to Bristol University for genotyping using the Axiom 35k SNP Wheat Breeders' Genotyping Array [66]. This array has previously been used for identifying QTLs and population structure in tetraploid wheat [67,68].
Marker allele calls were made using the Affymetrix Axiom Analysis Suite (www.thermofisher.com) with some alterations from the default parameters to account for the ploidy and diversity of the material: a dish quality control (DQC) threshold of 0.8, a genotype call rate threshold of 96% and a quality control (QC) call rate of 96% were used for the SNP calling. All markers with a call rate between 96% and 98% were manually inspected and recalled or disregarded appropriately. Markers named as 'no minor homozygote' were manually inspected and the missing homozygote classes were recalled where possible. Additionally, markers with a QC call rate between 95% and 96% were manually inspected and where clustering was still well defined, the markers were included in the analysis.
Further quality control and analysis was completed using R (version 3.5, [69]) in RStudio [70]. Markers that were monomorphic between the parents or heterozygous within a parental line were removed. Based on genotypes, Euclidean distances between individuals were calculated and plotted through principal coordinate decomposition, which was implemented by the R package ape [71]. Two pairs of individuals appeared to be more closely related than expected, suggesting possible out-crossing or mixed samples at some stage of the single seed decent program. One individual from each of these pairs was removed from the analysis. Heterozygosity frequency in markers and individuals were plotted and visually inspected, and a single outlier genotype and markers with over 25% heterozygote calls were removed. ChiSquare tests were used to remove markers that differed significantly from the expected allele frequency (1:1) typical of Mendelian inheritance. p-values (p) from these tests were adjusted for multiple tests using a Bonferroni correction. After the QC, 4274 markers and 112 individuals remained in the dataset.

Genetic Linkage Map Formation and QTL Mapping
The R package ASMap (V -1.0-4, [33]) was used for the map construction and remaining quality control. The 'pullCross' function was used to remove 2144 co-located markers. Using the 'mstmap' function, a numerical threshold value of p = 1 × 10 −12 was set for linkage group clustering, which separated the markers into an initial 25 different linkage groups. Using a one-tailed test of a Poisson mean and a Bonferroni adjustment for multiple testing, eight individuals were removed with a significantly different (p < 0.01) crossover number from the population median. Single markers that formed a unique linkage group or which had more than two double crossovers were removed from the map. Following Broman [72], genotype error rate was estimated through maximum likelihood as 0.0025 and genotypes with high error logarithm-of-odds (LOD) scores were identified and removed (LOD > 6). Markers from the SNP array were all run through the Basic Local Alignment Search Tool (BLAST+) command-line applications [73] against the RefSeq v1.0 wheat genome assembly [35] using default parameters. All markers with a single hit to the A genome were selected and retained for analysis of the A genome regardless of other genome hits; similarly all markers with a single hit to the B genome were selected for analysis of the B genome, regardless of other genome hits. The newly formed linkage groups were then compared to the physical positions of the extracted markers from the BLAST+ results and the breeders' array consensus map [66]: 28 linkage groups were merged and renamed appropriately. Marker physical positions were plotted against the mapped genetic position for each linkage group and 21 markers were dropped that were mapped to different physical and genetic chromosomes. Finally, heat maps for each chromosome were generated to validate marker ordering and check for linkage outside of the central diagonal. No issues were found through the heat maps and no changes were made. The completed map was formed from 2066 markers and 104 individuals.
QTL mapping was carried out for every trait using R/QTL (V -1.44-9, [34]) and the map consisting of 2066 markers. The genetic marker data used is shown in Supplementary Table S4. Phenotype data was used from 111 and 108 individuals from the 2017 field trial and the 2019 pot trial, respectively. To avoid violating the underlying assumptions of the model and creating false positives, a late flowering outlier was removed from the 2019 F T data before the QTL analysis was completed, the data point was a BLUE estimated from only a single replicate as the other replicate failed to grow. Following Broman and Sen [36], QTL genotype probabilities were calculated using the 'calc.genoprob' function with a step size of 1 cM. A two-stage QTL analysis was used. Firstly, a single genome scan was completed using the 'scanone' function. The LOD thresholds at significance levels of 5% and 10% (0.05 and 0.1 Alpha) were determined through running 5000 permutations within the 'scanone' function. If a QTL passed the determined 10% threshold, it was advanced to the second stage of analysis, where the functions 'makeqtl' and 'fitqtl' were used to fit a multiple-QTL model. For each trait, the function 'addqtl' was used to scan for additional QTLs, which were added to the model if the resulting LOD score from the function was over the initial 5% threshold. Where multiple QTLs were found, interactions between all QTLs in each model were tested for using the 'addint' function. QTL positioning was refined using the function 'refineQTL' and 1-LOD support intervals were determined using the function 'lodint' and were expanded to the nearest marker. Haley-Knott regression was implemented for all mapping functions. The physical positions on IWGSC Refseq v1.0 for all peak and 1-LOD support interval flanking markers were obtained from the previous BLAST+ results.

Conclusions
A novel T. dicoccum mapping population was created and QTLs were successfully identified using a map that we created consisting of 2066 SNPs. Candidate QTLs for A L and F T were mapped consistently over different trials and were likely linked to the major locus Hd and Vrn-B3. Novel candidate QTLs were potentially identified for FL W , T D and FL L but further work is needed for validation. Physiological trait compensations and the environment appeared to influence S D and no QTLs were found. The confidence intervals surrounding QTLs were relatively large and fine mapping is needed to provide improved accuracy on physical locations, which would also provide clarity on the putative relationship between the identified QTLs and major loci.
Supplementary Materials: The following are available online at http://www.mdpi.com/2223-7747/9/7/829/s1: Figure S1: Ear diversity within the PS1 population and photos of the mapping parents, Figure S2: Frequency histograms of the six traits across both years, Table S1: Summary statistics and two sample t-test results using the mapping parent trait data, Table S2: Summary statistics and two-sample Wilcoxon test results using the mapping parent trait data.