A population of 111 F10-derived pea recombinant inbred lines from the cross Dark Skin Perfection × 90–2131 was produced by single-seed descent at USDA-ARS, Pullman, WA, USA. DSP (Dark Skin Perfection, Unilever Limited 1960, also designated as W6 17516) is a spring pea cultivar used for freezing and canning, with white flowers, normal leaves, wrinkled seeds with white hilum and it is susceptible to A. euteiches. 90–2131  is a pea germplasm with white flowers, normal leaves, dimpled seeds with black hilum and it is partially resistant to A. euteiches in France and in the USA. The parental lines, 90–2131 and DSP, were used as check lines in field- and controlled-condition disease resistance tests of the DSP × 90–2131 RIL population. The three RIL parental lines previously studied in , Baccara (susceptible), PI180693 and 552 (partially resistant), were also used in the controlled condition assays. The MN313 line was included in the assays for distinguishing the two main pea pathotypes of A. euteiches. The pea spring variety Solara (susceptible) was used as the adjacent control in the field assays of the RIL population.
Six pure culture strains of A. euteiches virulent on pea were used for resistance evaluation of the DSP × 90–2131 RIL population in controlled conditions. These strains were RB84, Ae106, Ae78, Ae85, Ae87 (referred to as strain SP7 in ) and Ae109 (referred to as strain Ae467 in  and ). The six strains were chosen based on their various geographical origins and pathotype groups (Additional file 1). The strains Ae106 and Ae87, were previously used for disease screening of the Puget x 90–2079 RIL population  and the strains RB84 and Ae109 were used in the Baccara × PI180693 and Baccara x 552 RIL populations screenings .
Disease resistance evaluations
In controlled conditions, resistance of the DSP × 90–2131 RILs to pure-culture strains of A. euteiches was assessed on 14-day-old seedlings in a climatic chamber, as described in , with one strain and four blocks per test. Each block included all RILs and control lines, with five plants/line/block. A root rot index (RRI) ranging from 0 to 5 was calculated as the mean disease severity score on the five plants per line, as described in .
In the field, the DSP × 90–2131 RILs were evaluated over an international Aphanomyces infested nursery network, described in , over four years (2000 to 2003) and five locations in France (Riec-sur-Belon, Finistère (RI); Dijon-Epoisses, Côte d’Or (DI); Templeux, Somme (TPX)) and in the United States (Pullman, WA (PLM); LeSueur, MN (LS)) [in 2000: PLM, LS; in 2001: LS; in 2002 : RI, DI, LS; in 2003: RI, DI, TPX, PLM, LS]. Field assays were carried out using experimental designs as described in . For field experiments conducted in France (two years, three locations) and in Pullman in 2003, the design included a check plot of the susceptible cultivar Solara or DSP every two to four plots, in order to adjust disease severity scores for local disease variations in the soil  using the formula described in [46, 77]. Two disease criteria were used to assess resistance for each plot, as described in : (i) the root rot index (RRI), using a 0–5 scoring scale, evaluated each year in French nurseries, except at DI in 2002, and (ii) the aerial decline index (ADI) evaluated once, twice or three times in all USA and French disease nurseries (except at RI in 2002), using a 1–5 and 1–9 disease scoring scale, respectively.
Evaluation of morphological and phenological traits
The DSP × 90–2131, Baccara × 552 and/or Baccara × PI180693 RIL populations were evaluated for two agronomic traits: earliness at flowering (FLO) and plant height (HT).
The FLO trait was evaluated in the Aphanomyces-infested nursery of DI in 2003 (DSP × 90–2131 RILs), 2004 (Baccara × 552 RILs), 2006, 2007 and 2008 (Baccara × PI180693 RILs) using the experimental design established for resistance evaluation. The FLO trait was also evaluated in a healthy nursery at Rennes-Le Rheu (Ille-et-Vilaine, FR (REN)) in 2002 (DSP × 90–2131 RILs), 2005 and 2008 (Baccara x PI180693 RILs) using a randomized complete block design with 1 and 3 block(s) in 2002/2005 and 2008, respectively (40 plants/plot in a two m-long twin rows). The FLO trait was scored on each plot as the number of days to 50% bloom (FLO1) or to 100% bloom (FLO2) from the first day of the year. The HT trait was evaluated in the same REN healthy nurseries as used for FLO evaluation in 2002 (DSP × 90–2131 RILs), by measuring the average height of five plants at maturity in a whole plot.
Molecular markers and genetic mapping
The DSP × 90–2131 RIL population was genotyped using simple sequence repeat (SSR) from , random amplified polymorphic DNA (RAPD)  and known-function genes [48, 49]. Two morphological traits were also scored: Pl for hilum colour and R for round/wrinkled seeds. DNA extractions and PCR amplifications were performed as described in  and in [48, 49, 53]. The Puget × 90–2079 RIL population was genotyped using additional SSR, RAPD, SCAR and known-function gene markers from [49, 53, 69, 78–80] compared to markers reported in . Marker coding and the genetic map (in cM Kosambi) from the DSP × 90–2131 RIL population were established with a minimum LOD score threshold of 3.0 and a maximum recombination frequency of 0.4, as described in . For each locus, adjustment of allelic segregation to the expected 1:1 Mendelian ratio was analyzed using a χ2 test (α=0.01). Additional markers were placed on the framework Puget × 90–2079 genetic map reported in , using the “assign” and “place” commands of MAPMAKER/EXP version 3.0b .
Statistical and QTL analyses
Statistical and QTL analyses were conducted from field and controlled-condition disease scoring data obtained in the DSP × 90–2131 RIL population and from earliness and plant height data obtained in the DSP × 90–2131, Baccara × PI180693 and Baccara × 552 RIL populations, for each scoring variable in each environment. Statistical analysis of each data set was carried out using a two-way ANOVA estimating genotype and block effects and normality of residuals was analyzed, as described in . Broad sense heritability, RIL least-square means used for linkage analysis and Pearson correlation coefficients (r
) between adjusted mean data were also estimated as described in .
Additive-effect QTL analysis from each data set was performed by Composite Interval Mapping  using Windows QTL Cartographer 2.5 software , as described in . Using the permutation test with 1000 permutations, minimum LOD thresholds of 2.9 (for DSP × 90–2131 RIL population) and 2.8 (for the Baccara × PI180693 and Baccara × 552 RIL populations) were chosen for all the traits to declare a putative QTL significant, corresponding to a genome-wide α error risk of 5%. Two QTL were considered as belonging to the same genomic region when their one-LOD drop-off confidence intervals overlapped. QTL for resistance to A. euteiches, earliness and plant height were named “Ae-Ps”, “Flo-Ps” and “HT-Ps”, respectively, followed by the linkage group number and the QTL number within the linkage group for each trait. Based on common markers between genetic maps, Aphanomyces resistance QTL common to those previously published [31, 46] were named as described in . The most significant pairwise epistatic interactions were detected for each resistance variable between all possible marker pairs of the DSP × 90–2131 genetic map, as described in  using a detection threshold of P < 7,1.10-6 and R
2 > 5%.
From the Puget × 90–2079 RIL population, additive-effect QTL were re-detected, as described in , from the updated genetic map generated in this study and the phenotypic data reported in [31, 45].
Using all the additive-effect QTL identified for resistance to A. euteiches, earliness and plant height from Puget × 90–2079 [31, 78], Baccara × PI180693, Baccara x 552  and DSP × 90–2131 (this study) RIL populations, a QTL meta-analysis was performed using the MetaQTL software version 1.0 .
QTL meta-analysis was conducted in three steps, according to details given in . For each QTL, data given to the software were the QTL position (LG, position on the LG at the LOD peak), the upper and lower bound and LOD decrease of the QTL confidence interval, the percentage of variation (R
) individually explained by the QTL, the trait related to the QTL and the size of the mapping population used for the QTL detection. First, a single consensus marker map was built by integrating the available genetic maps, based on common markers designated with common names between the maps, using the ConsMap command of the software. The implemented method applied a weighted least-square strategy, using individual distances between markers in each individual map, to determine marker order and position on the consensus marker map. The InfoMap command was carried out in order to list markers whose orders were not consistent between the different individual maps. Markers with inconsistent positions were removed from the consensus marker map.
Second, the QTL detected in each study were projected onto the consensus map, using the QTLproj command of the software. This command enabled the homothetic projection of the individual QTL positions and confidence intervals based on a scaling rule between QTL-flanking marker positions on the individual maps and on the consensus map. QTL projection was carried out using LOD-1 confidence intervals of all individual QTL, for graphical representation and for identifying main genomic regions comprised of overlapping QTL intervals. For the meta-analysis, the projection was carried out using independent individual QTL, since the QTL meta-analysis algorithm implemented in the software assumes that the input mapping studies are independent from each other. Independent QTL were selected as follows. In a given mapping population and for a given variable, when QTL detected from various environments or strains had overlapping confidence intervals, only the QTL with the greatest proportion of the phenotypic variation (R
2) and the smallest confidence interval (if highest R
2 were equals) was retained. For example, among the five individual QTL detected at the top of LGII for the ADI variable in different environments (RI04, RI05, PLM04, TPX04) in the Baccara x 552 population, the one having the greatest R
(15%) was kept for the meta-analysis and the four others were removed. Confidence intervals (CI) of independent QTL were estimated from the QTL R
2 values, using the empirical formula proposed by :
, where N is the population size. This formula usually gives larger confidence intervals than the usual interval length of LOD-1 decrease.
Finally, the QTL meta-analysis algorithm, computed in the QTLClust command of the software, was used to determine the most likely number of meta-QTL on a given chromosome and to estimate their corresponding positions and confidence intervals. Meta-QTL considers QTL positions and corresponding confidence intervals in individual experiments, after projection onto the consensus map. We used the --cimode option 4 of QTLClust command that consider the maximum confidence interval between the LOD-1 decrease value reported and the QTL R
value estimated by the formula of . On a given chromosome, projected independent QTL were clustered into all possible numbers of hypothetic clusters of meta-QTL (K), for which a Gaussian mixture model estimates meta-QTL positions and confidence intervals . The optimal K was determined using the highest weight of evidence values estimated for five information-based criteria , which were computed for each K. The most frequent optimal K value given by the different criteria was selected, which always corresponded to the optimal K value determined by the Akaike Information Criterion. The position, 95% confidence interval and probability of individual QTL belonging to the meta-QTL, were given by the software for each meta-QTL.
Two meta-analysis runs were conducted separately for Aphanomyces resistance QTL data and earliness and plant height QTL data, respectively, resulting in the detection of meta-QTL for two types of traits, which were then compared.
Meta-QTL for resistance to A. euteiches and meta-QTL for earliness and plant height were named “MQTL-Ae” and “MQTL-Morpho”, respectively, followed by the meta-QTL number on the whole pea genome.
Identification of candidate genes underlying meta-QTL
Positional candidate genes included within the support interval of main resistance meta-QTL regions were identified using the Pea-M. truncatula toolkit of Bordat et al.. Mining the high colinearity between the pea and M. truncatula genomes, the toolkit allowed placing in silico 5460 pea Unigenes on the pea consensus map of , from positions of their best homologs on the M. truncatula genome. Positions of main meta-QTL regions identified in this study were estimated on the pea consensus map of  based on common markers, mainly SSRs. On the pea map of , intervals between adjacent marker genes covering main meta-QTL regions were identified and a list of Unigenes contained in these intervals was established. Intervals defined from adjacent marker genes genetically mapped in  were often larger than the reduced confidence interval of meta-QTL or absent. In these cases, intervals at most probable meta-QTL positions were chosen. Positional candidate genes were examined to identify those known to be involved in disease resistance in plants.