Mapping QTL associated with partial resistance to Aphanomyces root rot in pea (Pisum sativum L.) using a 13.2 K SNP array and SSR markers

A stable and major QTL, which mapped to an approximately 20.0 cM region on pea chromosome 4, was identified as the most consistent region conferring partial resistance to Aphanomyces euteiches. Aphanomyces root rot (ARR), caused by Aphanomyces euteiches Drechs., is a destructive soilborne disease of field pea (Pisum Sativum L.). No completely resistant pea germplasm is available, and current ARR management strategies rely on partial resistance and fungicidal seed treatments. In this study, an F8 recombinant inbred line population of 135 individuals from the cross ‘Reward’ (susceptible) × ‘00-2067’ (tolerant) was evaluated for reaction to ARR under greenhouse conditions with the A. euteiches isolate Ae-MDCR1 and over 2 years in a field nursery in Morden, Manitoba. Root rot severity, foliar weight, plant vigor and height were used as estimates of tolerance to ARR. Genotyping was conducted with a 13.2 K single-nucleotide polymorphism (SNP) array and 222 simple sequence repeat (SSR) markers. Statistical analyses of the phenotypic data indicated significant (P < 0.001) genotypic effects and significant G × E interactions (P < 0.05) in all experiments. After filtering, 3050 (23.1%) of the SNP and 30 (13.5%) of the SSR markers were retained for linkage analysis, which distributed 2999 (2978 SNP + 21 SSR) of the markers onto nine linkage groups representing the seven chromosomes of pea. Mapping of quantitative trait loci (QTL) identified 8 major-effect (R2 > 20%), 13 moderate-effect (10% < R2 < 20%) effect and 6 minor-effect (R2 < 10%) QTL. A genomic region on chromosome 4, delimited by the SNP markers PsCam037549_22628_1642 and PsCam026054_14999_2864, was identified as the most consistent region responsible for partial resistance to A. euteiches isolate Ae-MDCR1. Other genomic regions important for resistance were of the order chromosome 5, 6 and 7.


Introduction
Field pea or "dry pea" (Pisum sativum L.) is an economically important cool-season legume crop that is cultivated widely in different parts of the world (Hossain et al. 2012).
Unfortunately, the production of field pea is affected adversely by the pea root rot complex (PRRC) (Bailey et al. 2003;Xue 2003;Chang et al. 2013;Wu et al. 2018a, b). The soilborne oomycete Aphanomyces euteiches Drechs. is a dominant pathogen in the PRRC (Jones and Drechsler 1925). This pathogen is favored by saturated soil conditions and poor drainage. The oospores can survive in the soil for up to 10 years (Papavizas and Ayers 1974;Holliday 1980). Under conducive environmental conditions, yield losses in pea as high as up to 86% can occur in fields infested with A. euteiches (Pfender and Hagedorn 1983). In Canada, Aphanomyces root rot (ARR) outbreaks have been reported only Communicated by Janila Pasupuleti.
1 3 recently, either because pea production in the same fields over multiple years resulted in a buildup of the pathogen, or because symptoms of ARR can now be better distinguished from other root rot pathogens (Hwang and Chang 1989;Xue 2003;Chatterton et al. 2015;Wu et al. 2017). Aphanomyces root rot is characterized by the formation of soft and water-soaked rootlets with a honey-brown or blackish-brown color. Reductions in seedling emergence and seedling blight have also been reported. As infected plants continue to grow, secondary infection causes development of brown lesions and cortical decay of the belowground tissues. The uptake of water and nutrients in diseased plants also is reduced, which can result in wilting and death of the plants (Chatterton et al. 2015).
Chemical strategies appear to be of limited value in the control of ARR, due to a lack of effective commercial fungicides (Pilet-Nayel et al. 2002;Wu et al. 2019). The fungicides Apron FL (metalaxyl-M) and fosetyl-AL were reported to increase seedling emergence and delay the infection of field pea by A. euteiches in the early growth stages, but they failed to suppress ARR in the later stages, leading to yield losses (Oyarzun et al. 1990;Xue 2003). The seed treatment INTEGO™ Solo, containing the chemical ethaboxam (Valent Canada, Inc. Guelph, Ontario), is the only product registered in Canada found to suppress the growth of A. euteiches. Chemical fungicides, however, are host non-specific and hence may affect beneficial microorganisms in the soil as well as pollinators. Furthermore, fungicides may leave residues on crops and persist in the environment. Cultural disease management methods, such as longer rotations with non-host crops and the avoidance of infested fields, have had some success but are not always practical (Malvick et al. 1994;Conner et al. 2013). Genetic resistance could be the most promising approach to control ARR. However, pea cultivars completely resistant to ARR are not available (Pfender et al. 2001;Conner et al. 2013;Gossen et al. 2016). As a result, genotypes with only partial polygenic resistance are being used for the economic and durable control of ARR (Palloix et al. 2009;Kou and Wang 2010;Desgroux et al. 2016;Lavaud et al. 2015). Shehata et al. (1983) identified tolerance to ARR in some plant introduction (PI) lines of pea.
Previous studies to map quantitative trait loci (QTL) utilized morphological traits (e.g., leaf morphogenesis, hilum color on seeds and anthocyanin production), isozymes, and amplified fragment length polymorphism (AFLP), random amplified polymorphic DNA (RAPD), simple sequence repeat (SSR), inter-simple sequence repeat (ISSR) and sequence-tagged site (STS) markers (Pilet-Nayel et al. 2002Hamon et al. 2011Hamon et al. , 2013Lavaud et al. 2015). The number of PCR-based markers used for genotyping in these studies has ranged from 150 to 350. As such, the confidence interval for the detected QTL is often very large and the identified QTL usually contain a small number of markers. The use of high-density SNP arrays or next-generation sequencing technology (NGS) in peas has only been reported recently (Sindhu et al. 2014;Tayeh et al. 2015;Desgroux et al. 2016;Gali et al. 2018). Reducing the confidence interval in the QTL regions will be required to identify associated markers for marker-assisted selection (MAS). The availability of the pea reference genome (Kreplak et al. 2019) will facilitate the development of different marker types for genotyping this important crop. Pilet-Nayel et al. (2002 identified one major QTL Aph1 and three minor QTL Aph9, Aph10 and Aph11 located on pea chromosome IV to be associated with resistance to ARR in a recombinant inbred line (RIL) population. Additionally, eight minor QTL, Aph2 located on chromosome V, Aph3, Aph4 and Aph5 on chromosome I, Aph8 on chromosome III, Aph12 on chromosome VI, and Aph6 and Aph13 on chromosome VII, were found to be associated with the resistance to ARR in the same RIL population (Pilet-Nayel et al. (2002. The major QTL (Aph1) accounted for up to 47% of the variability in partial resistance. Hamon et al. (2011Hamon et al. ( , 2013 identified five highly stable QTL associated with partial resistance to ARR in pea. The QTL Ae-Ps7.6 located on chromosome VII had a major effect on resistance and explained up to 42.2% of the phenotypic variation in 32 of 37 disease variables in two RIL populations. Four other QTL, namely Ae-Ps1.2 on chromosome I, Ae-Ps2.2 on chromosome II, Ae-Ps3.1 on chromosome III, and Ae-Ps4.1 on chromosome IV, accounted for up to 14.4, 26.9, 29.9 and 24.5% of the phenotypic variation, respectively, in 13, 22, 11 and 14 of 37 disease variables in the same two RIL populations. Lavaud et al. (2015) identified Ae-Ps7.6 and Ae-Ps4.5, located on chromosomes VII and IV, respectively, as major-effect QTL for tolerance to ARR in near isogenic lines (NIL) of pea with different genetic backgrounds. The QTL Ae-Ps5.1, located on chromosome V, contributed the least to the resistance in the NILs (Lavaud et al. 2015). Desgroux et al. (2016) identified 52 QTL (with small-size intervals on all seven chromosomes) in a genome-wide association study (GWAS) of 175 pea lines. Collectively, these studies suggest that the major QTL for tolerance to ARR in pea are located on chromosomes IV and VII, with several minor QTL on chromosomes I, II, III and V (Pilet-Nayel et al. 2002Hamon et al. 2011Hamon et al. , 2013Lavaud et al. 2015;Desgroux et al. 2016).
The purpose of this study was to identify QTL for partial resistance to ARR using high-density SNP markers and SSR anchor markers in an F 8 RIL population of pea derived from the cross 'Reward' (susceptible) × '00-2067' (tolerant). In addition, QTL for two disease-related traits (foliar weight and vigor) and one agronomic trait (plant height) were evaluated. Lastly, the stability of the genetic loci controlling these four traits was examined under field and greenhouse conditions.

Fungal isolate
The A. euteiches isolate Ae-MRDC1, obtained from soil collected from an AAR disease nursery at the Agriculture and Agri-Food Canada (AAFC) Morden Research and Development Centre (MRDC), Morden, MB (lat. 49°11' N, long. 98°5' W), was used in the greenhouse experiments. This isolate was identified as virulence type I following . In brief, the oospore inoculum was generated in an oat broth as described by Papavizas and Ayers (1974). A preliminary inoculum density experiment was conducted to determine the concentration of inoculum that produced a differential reaction in the parents. This was done by first estimating the oospore concentration with a hemocytometer and adjusting the final concentration to 1 × 10 2 , 1 × 10 3 , 1 × 10 4 , 1 × 10 5 , and 1 × 10 6 oospores mL −1 with sterile deionized water. Based on this preliminary trial, an oospore concentration of 1 × 10 5 oospores mL −1 was selected for use in screening the parents as well as the RILs. The isolate Ae-MRDC1 was re-isolated from infected root tips and was confirmed to be A. euteiches by PCR analysis with the species-specific primer set, 136F/211R (Vandemark et al. 2002).

Phenotyping under field conditions
Field evaluation of the RILs was conducted in the naturally infested ARR disease nursery at the AAFC MRDC in 2015 and 2016. Pea genotypes planted in this nursery, which is situated on a loamy clay soil, have consistently developed severe ARR symptoms over many years, confirming strong disease pressure (Corner et al. 2013). Field layouts of the 135 RIL and three repeats of the two parental checks, '00-2067' and 'Reward', were generated as generalized lattice designs with the experimental design software CycDesigN® (VSN International 2015). The layout differed slightly between 2015 and 2016. In 2015, each replicate consisted of nine rows by 16 plots, with a check in each row and the two checks occurring once in each set of three rows. In 2016, the lattice layout was Latinized to account for any gradients from left to right, and up and down the field. The three replicates each with six rows were stacked up the field, and each row contained 24 plots formed by three blocks of eight plots (i.e., 6 rows × 3 blocks × 8 plots/replication). Two blocks each with 48 plots from each of the three replicates created three super-blocks. Each super-block contained 135 RILs and 3 repeats of the two checks.
Plots with 15 seeds/row and 30 seeds/row were sown on May 7, 2015, and May 9, 2016, respectively. Fertilizer applications and weed control were based on standard recommendations for field pea production in the region (Saskatchewan Pulse Growers 2000). Ten plants were uprooted from each plot for trait measurements. Root rot disease severity (DSF) and vigor (VF) were measured when the plants were at the pods setting stage on July 22, 2015, and July 19, 2016. Disease assessment (DSF) was carried out on a 0-9 scale following Conner et al. (2013), where 0 = healthy root system and 9 = roots totally destroyed and rotten. Vigor (VF) was scored on a 0-4 scale, where 4 = completely healthy plant with no wilting and 0 = stunted plants with yellow leaves and dead stem and leaves (Conner et al. 2013). The foliar weight (Fwt) of the plants evaluated for DSF and VF was obtained after oven drying at 30 °C for 10 days.

Phenotyping under greenhouse conditions
Eight seeds of each RIL and both parents were germinated on moistened Whatman No. 1 filter paper in Petri dishes, with the seedlings transplanted into 7 × 7 × 10 cm plastic pots (1 plant/pot) filled with a sterilized potting mixture (Cell-TechTM, Monsanto, Winnipeg, MB) after 5 days. Briefly, a 5-cm-deep hole was made in the soiless mix in each pot, with one seedling placed in the hole and inoculated with a 1-mL aliquot of the oospore suspension prepared as described above. Non-inoculated plants of the parental cultivars 'Reward' and '00-2067' were included as controls. The rootlets of the inoculated seedlings were then covered with the potting mix, and the plants were placed in a greenhouse with a 12-h photoperiod at day and night temperatures of 22-28 °C and 15-18 °C, respectively. The pots were arranged in a randomized complete block design (RCBD). Plants were watered daily in the morning and evening to ensure high moisture levels in the 1 3 potting mixture. The entire inoculation experiment was run independently three times, referred to as greenhouse experiments 1 (GH1), 2 (GH2) and 3 (GH3).
Plant height (HGH) was measured from the top leaf to the soil line at the end of the second week after inoculation. At four weeks, the plants were carefully uprooted, washed under standing water, and assessed for disease severity (DSGH), plant vigor (VGH) and dried foliage weight (DFGH) as described for the field study. All parameters were measured by each replicate (each single plant).

Statistical analysis of phenotypic data
All of the variables from the field and greenhouse trials were analyzed by station-year or single greenhouse experiment and pooled conditions (i.e., year × location) using R software (R core team 2019). The model for single environment analysis was: Pij = μ + Gi + Rj + eij, where Pij is the score of the ith RIL located in the jth replicate, μ is the mean of all the data in a single site-year or greenhouse experiment, Gi is the ith RIL effect or genetic effect, Rj is the jth replicate effect or blocking effect, and eij is the residual variance. The model for multiple station-years or greenhouse experiments was: Pijk = μ + Gi + Rj + Lk + GLijk + eijk, where Pij was the score of the ith RIL located in the jth replicate, μ the mean of all pooled data, Gi the ith RIL effect or genetic effect, Rj the jth replicate effect or blocking effect, Lk the kth environment or location effect, GLijk the interaction of RIL effect and environment, and eijk was the residual variance. The best linear unbiased predictors (BLUPs) were also calculated using the package 'Phenotype' in R software (R core team 2019), since the year-station effects were significant. The entry-based broad-sense heritability (H 2 = σ 2 G /σ 2 P ) was calculated as h 2 = σ 2 G /[σ 2 G + (σ 2 e /r)] for single site-years or greenhouse experiments for which σ 2 G is the genetic variance, σ 2 e the residual variance, k the number of environments, and r the number of replicates per line. The heritability of BLUPs was also obtained using package 'Phenotype' in R software (R core team 2019). The Pearson correlation coefficient of the means in each environment as well as BLUPs was calculated for each variable in the different single site-years or greenhouse experiments. Correlation analysis of the Pearson correlation coefficient among variables also was used to evaluate the relationship between disease tolerance and the agronomic variables. The Shapiro-Wilk test was used to test for normality in the phenotypic data. The power to detect a QTL for RILs of 135 individuals was determined using 'powercalc' in R/qtl for all the traits using the total data from the field and greenhouse, independently. The biological replicates, genetic variance and environmental variance were extracted and calculated from ANOVA.

Genotyping with SNP and SSR markers
The 135 RILs and the parents were first genotyped at Trait-Genetics GmbH, Gatersleben, Germany, using 13,204 highquality SNP markers selected from a 248,617 SNP marker set developed by Tayeh et al. (2015). The 13,204 SNPs were all derived from gene-encoding sequences and were well distributed across the pea genome (Tayeh et al. 2015). Filtering was conducted to remove failed SNP reactions, markers lacking polymorphism in the parents, low coverage site markers, markers with MAF ≤ 0.05, markers missing data for > 5% of the accessions, and those with segregation distortion.
Two hundred and twenty-two microsatellite markers, reported by Loridon et al. (2005) to be well distributed along the seven linkage groups of pea, were also used to genotype the 135 RILs and the parents. PCR assays were carried out in a 12-µL reaction mixture containing 20 ng of genomic DNA, 1 × Taq buffer, 2.0 mM MgCl 2 , 200 µM dNTPs, 0.4 µM forward primer modified at the 5'-end with an M13 tail, 0.4 µM reverse primer, 0.2 µM fluorescently labeled M13 primer and 1.25 U Taq polymerase (Promega, Madison, USA). Amplifications were carried out in a Mycycler Thermal Cycler (Bio-Rad, Mississauga, ON, Canada) with 35 cycles of denaturation at 94 °C for 30 s (5 min for the first cycle), annealing for 45 s at a temperature based on the primers used, and extension at 72 °C for 1 min.
An aliquot of the PCR products was separated by capillary electrophoresis on an ABI PRISM 3730xl DNA analyzer (Applied Biosystems, Foster city, CA). In addition, the amplified products of the polymorphic markers were separated by electrophoresis on 8% polyacrylamide gels (PAGE) at 150 V for 2 h. The amplified fragments were stained with silver nitrate and photographed with a UV Transilluminator (Bio-Rad Laboratories, Mississauga, ON, Canada).

Linkage map construction
A linkage map was constructed with the filtered SNP and the polymorphic SSR markers following the two-step mapping strategy of Perez-Lara et al. (2016). In brief, a draft linkage map was generated using minimum spanning tree map (MSTMap) software (Wu et al. 2008) with a strict cutoff p value of 1E −10 and a maximum distance between markers of 15 cM. Only one of multiple markers in the same positions was retained for linkage analysis. The draft map was then refined using MAPMAKER/EXP 3.0 (Lincoln et al. 1992) with a logarithm of odds (LOD) score ≥ 3.0 and recombination fraction (θ) value ≤ 0.40. The Kosambi map function (Kosambi 1944) was used to calculate the genetic distances (in cM) between the markers. The linkage groups were assigned to chromosomes based on the consensus SNP map of pea developed by Tayeh et al. (2015). Genetic linkage maps were constructed with MapChart v. 2.32 (Voorrips, 2002).
Additive-effect QTL were detected by Composite Interval Mapping (CIM) using WinQTL Cartographer v2.5 with the following parameters: 1 cM walking speed, forward and backward regression method, window size 10 cM, five background cofactors, 1000 permutations and P < 0.05 (Wang et al. 2012). The QTL positions were identified at regions where the LOD score reached values ≥ 3.0. The confidence interval of each QTL was defined by the consensus region bordered by the five environments (two field and three greenhouse experiments).
Additive-effect QTL were named according to Tayeh et al. (2015) except that the name of the A. euteiches isolate (Ae-MRDC1) was indicated, to avoid confusion, for example, of Ae MRCD1 -Ps2.1 (this study) with Ae-Ps2.1 (Tayeh et al. 2015). A similar nomenclature was adopted for foliar weight (Fwt-Ps2.1), vigor (Vig-Ps2.1) and plant height (Hgt-Ps2.1), where "Ps" = Pisum sativum, the first number = pea linkage group (Tayeh et al. 2015), and the second number = the serial number of the QTL on the linkage group. The location of QTL on the pea chromosome and pseudomolecules were in accordance with Neumann et al. (2002) and Kreplak et al. (2019), respectively. The QTL were classified as stable if they were confirmed in at least two of five environments. The percentage of phenotypic variation explained (PVE) due to a particular QTL was estimated by the coefficient of determination (R 2 ). Furthermore, QTL with R 2 > 20%, 10-20% and < 10% were arbitrarily classified as major, moderate or minor-effect QTL, respectively.
The origins of favorable alleles for individual traits were assigned to different parents following Lubberstedt et al. (1997) and Zaidi et al. (2015). Alleles coming from '00-2067' were coded as "2", while alleles from 'Reward' were coded as '0". The additive effects of each QTL were calculated by deducting the phenotypic average of all individuals carrying the "0" allele from that of individuals with the "2" allele. A negative sign of the additive effect for root rot severity indicated that the favorable allele for the traits originated from the parent '00-2067', while a positive sign indicated it originated from 'Reward'. In contrast, a positive sign of the additive effect for foliar weight, vigor and plant height indicated that the favorable allele for these traits originated from '00-2067', while a negative sign indicated that the favorable allele originated from 'Reward'.

Epistatic-effect QTL analysis
Epistatic-effect QTL (QTL × QTL) were detected with Ici-Mapping V.4.1 using the ICIM-EPI method (Meng et al. 2015). The mapping parameters were the same as those used for the CIM above. Epistatic-effect QTL were named with the prefix "E" followed by the QTL name and a serial number (e.g., E.Ae MRCD1 -Ps1, E.Fwt-Ps1, E.Vig-Ps1 and E.Hgt-Ps1). The R 2 significance thresholds selected to classify epistatic-effect QTLs as major, moderate or minor were arbitrarily set at R 2 > 15%, 7.5-15% and < 7.5%, respectively.

Candidate gene prediction
The sequences of the SNP markers found to be associated with Aphanomyces root rot tolerance loci were used in BlastN searches of the Pulse Crop Database (www. pulse db. org/) to determine their possible functions. An E-value ≤ E-20 and a percentage identity of ≥ 95% were used as the thresholds to confirm the putative functions of candidate genes.

Phenotypic trait analysis
The evaluation and frequency distribution of the 135 RILs for disease severity, foliar weight, vigor and plant height for each of the two field seasons and three greenhouse experiments, as well as for the best linear unbiased predictors (BLUPs) of the experiments, are presented in Table 1 and Fig. 1, respectively. The power to detect QTL given the 135 RIL population size used in this study was 0.99-1.00 for disease severity, vigor and foliar weight, and 0.67 for plant height.

Disease severity variation in the RIL population
The parental cultivar '00-2067' was tolerant to A. euteiches isolate Ae-MRDC1, with an estimated mean DSI of 2.3 ± 1.0 SE, 1.5 ± 1.1 SE and 1.7 ± 0.6 SE for the three greenhouse experiments (DSGH1, DSGH2 and DSGH3), respectively. The same cultivar was found to be moderately susceptible in the 2015 (DSF15) and 2016 (DSF16) field trials at Morden, MB, with estimated means of 6.1 ± 0.4 SE and 5.9 ± 0.6 SE, respectively. In contrast, the parent 'Reward' Table 1 Statistical summary of the traits for the parents (pea cultivars '00-2067' and 'Reward'), the RIL population and the Shapiro-Wilk test based on three greenhouse experiments, field experiments conducted in Morden, MB, in 2015 and 2016, and the best linear unbiased predictors (BLUPs) of the greenhouse and the field experiments *GH 1 = greenhouse experiment 1, GH 2 = greenhouse experiment 2 and GH 3 = greenhouse experiment 3 Trait was susceptible to ARR, with estimated means of 5.0 ± 1.3 SE, 6.1 ± 0.9 SE and 5.1 ± 0.5 SE for DSGH1, DSGH2 and DSGH3, and 8.1 ± 0.5 SE and 7.6 ± 1.1 SE for DSF15 and DSF16. ANOVA of each field and greenhouse (GH) test indicated a significant effect of genotype on the RIL population (P < 0.05). The year-station effect for root rot severity in Fig. 1 Correlation analysis of a all the single environment means and BLUPs of Aphanomyces root rot disease severity, pea dry foliar weight, vigor and plant height under field and greenhouse conditions, indicating the coincidence among single environments and BLUPs; and b relationship among all the traits in each single environment and BLUPs. The frequency distributions are shown in the bar graphs across the diagonal. The correlation coefficients and scatter plots between pairs are indicated above and below the diagonal, respectively. The significance levels are noted with asterisks, where *indicates P < 0.05; **indicates P < 0.01; and ***indicates P < 0.001 the greenhouse conditions and the disease nurseries was significant (P < 0.05) (Supplementary Table S1). The correlation among individual experiments and BLUPs in the greenhouse and field was significant (r = 0.36-0.88, P < 0.001, Fig. 1a). Therefore, the BLUPS were applied instead of the pooled data in the multiple model. The BLUPs of '00-2067' were calculated as 0.2 for the greenhouse data (DSGHB) and as 5.0 for the BLUPs field (DSFB) data. In the case of 'Reward', the BLUPs were 5.9 in the greenhouse and 8.2 in the field. Based on a t-test analysis, the differences in disease severity between the two parental genotypes were significant in both the greenhouse and field experiments (P < 0.05). The frequency distributions of the estimated means of the disease severity in the greenhouse and field were continuous (Fig. 1a), with low values for skewness and kurtosis (Table 1). The Shapiro-Wilk test of the disease severity data for normality was significant for all of the field and greenhouse experiments except for DSGH2. The genotypic effects of disease severity were significant (P < 0.05) in the greenhouse and field for the individual experiments. Entry mean-based heritability of disease severity was high, ranging from 60 to 92% (Supplementary Table S2), which indicated a strong genetic effect of tolerance to A. euteiches that was transmitted from '00-2067' to individuals in RIL population.

Foliar weight variation in the RIL population
Foliar weight for the parental cultivar '00-2067' had an estimated mean of 3.0 ± 0.7 SE, 2.6 ± 0.5 SE and 2.8 ± 0.5 SE in the three greenhouse experiments (DFGH1, DFGH2 and DFGH3), respectively, and 15.9 ± 11.7 SE and 5.6 ± 2.1 SE in the 2015 (DFF15) and 2016 (DFF16) field trials. In the case of 'Reward', the estimated means were 2.3 ± 0.8 SE, 1.8 ± 0.8 SE and 2.4 ± 0.5 SE for DFGH1, DFGH2 and DFGH3, respectively, and 2.7 ± 1.2 SE and 1.1 ± 0.5 SE for DFF15 and DFF16. The BLUPs for '00-2067' and 'Reward' were 2.8 and 1.8 for the greenhouse experiments (DFGHB), and 11.3 and 0.6 for the field trials (DFFB), respectively. Significant differences between the parental cultivars were found with respect to foliar weight in the field and greenhouse experiments (P < 0.05). In addition, significant RIL genotype effects were found in the ANOVA of foliar weight in the greenhouse (P < 0.05) but not in the field experiment. Foliar weight of the RIL population in both the greenhouse and field trials had a continuous frequency distribution (Fig. 1a). Based on the Shapiro-Wilk test, however, the data did not follow a normal distribution except for DFGH1 (Table 1). Foliar weight in the replicated greenhouse and field experiments was correlated significantly (r = 0.50-0.86, P < 0.01). Heritability of foliar weight was high, ranging from 67 to 94% (Supplementary Table S2). This suggested that a significantly high percentage of genotypic effect (P < 0.05) was transmitted from the parents to the RIL population.

Plant height variation in the RIL population
Differences in plant height between the parents '00-2067' and 'Reward' were significant (P < 0.05) in the greenhouse experiment 1 (HGH1). The estimated means in plant heights for '00-2067' were 223.6 ± 13.8 SE, 173.0 ± 25.6 SE and 164.2 ± 30.5 SE for HGH1, HGH2 and HGH3, respectively, while the estimated means for 'Reward' were 179.5 ± 28.7 SE, 140.3 ± 36.4 SE and 157.6 ± 23.1 SE for HGH1, HGH2 and HGH3, respectively. The BLUPs for plant height (HGHB) for '00-2067' were 183.7 mm, and for 'Reward' it was 130.1 mm. Significant RIL genotype effects were detected by ANOVA (P < 0.05) in all three greenhouse experiments. Frequency distribution of height in each greenhouse experiment was not normal based on the Shapiro-Wilk test ( Table 1). The height of both of the parents '00-2067' and 'Reward' was lower than the mean of the RIL population. The heritability of height was 93, 91, 92 and 94% for HGH1, HGH2, HGH3 and HGHB, respectively. (Supplementary Table S2).

Correlation analysis among traits
All of the traits in different environments were significantly and positively correlated with each other (0.33 < r < 0.95, P < 0.001). Correlation analysis among variables indicated 1 3 that all the traits were significantly correlated with each other in the individual experiments and in the BLUPs (0.31 < r < 0.90, P < 0.001), except for plant height, which was only coincidently correlated with vigor (0.18 < r < 0.38, P < 0.05) and dry foliar weight (0.59 < r < 0.71, P < 0.001). Root rot severity was negatively correlated with dry foliar weight and plant vigor in both the greenhouse and field experiments. This suggested that ARR had an adverse effect on plant growth (Fig. 1b).

Genetic linkage mapping
Filtering removed 10,154 (76.9%) of the SNP markers and 192 (86.5%) of the SSR markers used to genotype the 135 RIL population. Therefore, 3050 (23.1%) of the SNP and 30 (13.5%) of the SSR markers used for genotyping were retained for linkage analysis. The linkage analysis distributed 2999 (2978 SNP + 21 SSR) markers on nine linkage groups, while the remaining 81 (72 SNP + 9 SSR) markers were unlinked. The nine linkage groups represented all seven chromosomes of the pea genome ( Table 2). The length of the nine linkage groups ranged from 21.1 cM (linkage group 2) to 395.7 cM (linkage group 4). Linkage groups 1 and 2 from this study corresponded to linkage group I (Ia and Ib) of the pea genome by Tayeh et al. (2015) and to chromosome 2 by Neumann et al. (2002) ( Table 2). Linkage groups 3, 4, 5, 8 and 9 from this study corresponded to the linkage groups II, III, IV, VI, and VII, respectively, reported by Tayeh et al. (2015) and to chromosomes 6, 5, 4, 1 and 7, respectively, reported by Neumann et al. (2002) (Table 2). Linkage groups 6 and 7 from this study corresponded to linkage group V and chromosome 3 of the pea genome reported by Tayeh et al. (2015) and Neumann et al. (2002), respectively ( Table 2). The relationships between the linkage groups identified in this study and the pea reference genome reported by Kreplak et al. (2019) are provided in Table 2. The number of markers per chromosome in this study ranged from 103 in linkage group 8 (chrom 1/LG VI of the pea genome) to 334 in linkage group 9 (chrom 7/LG VII of the pea genome). The length of the seven chromosomes ranged from 155.6 cM in linkage group 8 (chrom 1/LG VI of the pea genome) to 395.7 cM in linkage group 4 (chrom 5/ LG III of the pea genome) and spanned a total length of 1704.1 cM. The marker density per cM ranged from 1.1 to 2.4 and averaged 1.8 markers per cM (Table 2). For convenience and ease of comparison with previous studies, we will use a modified linkage group designation of Tayeh et al. (2015) and the pea chromosome nomenclature of Neumann et al. (2002) in the rest of this paper.

Additive-effect QTL analysis
One thousand five hundred and seventy-seven (1577) of the 2999 markers (Supplementary Table S3) mapped to the same position as other markers and were excluded from the QTL analysis. Therefore, only 1422 (10.5%) of the initial 13,426 (13,204 SNP + 222 SS) markers used to genotype the 135 RIL population were used for QTL analysis. Collectively, 27 QTL related to tolerance to ARR, foliar weight, vigor and plant height were detected by the CIM method using WinQTL Cartographer v2.5 (Wang et al. 2012). Based on the R 2 values, 8 of these QTL were major-effect QTL, 13 were moderate-effect QTL and 6 were minor-effect QTL. Ten of the 27 QTL were detected Table 2 The distribution of single nucleotide polymorphism (SNP) and simple sequence repeat (SSR) markers on nine linkage groups representing all seven chromosomes of F 8 -derived recombinant inbred lines of the cross between the pea lines '00-2067' × 'Reward' α Pea chromosomes named according to Neumann et al. (2002), β Pea linkage groups named according to Tayeh et al. (2015) and γ pseudomolecule labels in the pea genome assembly v1a named according to Kreplak et al. (2019)  consistently in the greenhouse and field experiments, and could be classified as stable, while the remaining 17 were year-or experiment-specific QTL detected in only one environment (BLUPs data not counted).
The additive effects for Ae MRCD1

QTL for foliar weight
Seven putative QTL were detected by the CIM for foliar weight in the RIL population inoculated with A. euteiches isolate Ae-MRDC1 (Table 3). Fwt-Ps2.1 detected in the 2016 field (DFF16) experiment and the BLUPs of the field data (DFFB) data explained 5.3% and 6.5%, respectively, of the total variation in foliar weight.  (Fig. 3a).    (Table 3).

Epistatic interactions for QTL pairs
Three hundred and seventy putative digenic epistatic pairs were detected for root rot severity, dry foliar weight, vigor and plant height in the three greenhouse experiments and the two field experiments conducted in 2015 and 2016 1 3 (data not shown). The number of putative digenic interactions detected in the five environments ranged from 14 to 20, 17-27, 17-22 and 18-24 for root rot severity, dry foliar weight, vigor and plant height, respectively. Of the 370 putative digenic interactions, one of the QTL pairs had a PVE ≥ 15%, 19 had 7.5% ≤ PVE ≤ 15% and 350 had LG II (chrom 6); the moderate-effect QTL Vig-Ps2.1 was detected in the 2015 field experiment, while the minor effect QTL Vig-Ps2.2 was detected in the 2016 field experiment and the BLUPs of the field data. c Two major-effect QTL were detected on LG III (chrom 5); Vig-Ps3.1 was detected in the 2016 field experiment and in the greenhouse experiments 1 and 2, as well as in the BLUPs of the greenhouse data, while Vig-Ps3.2 was detected in the BLUPs of the field data. d The most stable major-effect QTL Vig-Ps4.1 and Vig-Ps4.2, which are located in close proximity on LG IV (chrom 4), were detected in the 2015 field experiment, in the greenhouse experiments 2 and 3, as well as in the BLUPs of the field data a PVE ≤ 7.5%. In the case of the BLUPs data, the number of putative digenic interactions detected was 20, 24, 25, 18, 21, 29 and 28 for DSFB, DSGHB, DFFB, DFGHB, VFB, VGHB and HGHB, respectively (Fig. 6). Of the 165 digenic interactions detected in the BLUPs analysis, 7 QTL pairs had 7.5% ≤ PVE ≤ 15% and 158 had a PVE ≤ 7.5%. A list of the single significant (PVE ≥ 15%) and 30 moderate (PVE ≥ 7.5%) digenic interactions is presented in Table 4.
Markers for 24 of the 27 major-and moderate-effect epistatic QTL (Table 4) were also linked to 17 additive-effect (moderate-effect QTL), Hgt-Ps4.1 (moderate-effect QTL), Hgt-Ps5.1 (minor-effect QTL) and Hgt-Ps7.1 (moderate-major-effect QTL) located on LG III (chrom 5), IV (chrom 4), V (chrom 3) and VII (chrom 7) were detected consistently in the greenhouse experiments and in the BLUPs of the greenhouse data Fig. 6 Epistatic QTL conferring tolerance to Aphanomyces root rot of pea in the pooled a greenhouse and b field data; dry foliar weight in the pooled c greenhouse and d field data; vigor in the pooled e greenhouse and f field data; and plant height in g the pooled greenhouse data, as indicated by QTL IciMapping V4.1 software. The dashed lines represent epistatic interaction pairs of epistatic QTL, while the numbers represent the LOD scores. α LGA-Linkage analysis according to the present study; β LG-Pea linkage groups named according to Tayeh et al. (2015) and γ Chrom-Pea chromosomes named according to Neumann et al. (2002)  These proteins are associated with plant defense responses and regulation, cell wall components and properties, biological processes, as well as to the response to abiotic and biotic stresses (Table S4).

Discussion
This study evaluated tolerance or partial resistance to ARR in an F 8 RIL population, with this resistance derived from the partially resistant parent '00-2067' (2013). Significant genetic effects within the RIL populations were detected for root rot severity, foliar weight, vigor and plant height. This was probably due to genetic differences in the parents (Conner et al. 2013), which were manifested as diversity alleles in the RIL population. The frequency distribution of the single experiments and the BLUPs data for all traits in the field trials or greenhouse experiments were continuous, but deviated from normality with various levels of skewness and kurtosis, which is not unusual for field disease data (Eskridge 1995;Feng et al. 2011;Coyne et al. 2015). This probably reflected environmental effects and the contribution of many different QTL, each of which was responsible for small increments in the resistance. Mohan et al. (1997) and Collard et al. (2005) reported that for a preliminary mapping study, a population size of 50-250 individuals was sufficient to reduce genotyping costs. Therefore, it was essential to determine the effective population size to obtain enough power for this analysis. Power analysis confirmed the sufficiency (close to 1) of the 135 RILs used in this study for disease severity, dry foliar weight and vigor. In contrast, the power of plant height (0.67) in the greenhouse was lower than the 0.8 threshold (Hu and Xu 2018;Kim 2016;Serdar et al. 2021). The small difference between parents for plant height indicated that this trait was not affected by ARR, but rather reflected the genetics of the field pea, which was also indicated by correlation analysis. The number of individuals included this study was within the range of 111-178 used by Pilet-Nayel et al. (2002 and Hamon et al. (2011Hamon et al. ( , 2013 for the detection of QTL associated with ARR resistance in field pea. Transgressive segregation, in which some lines were more resistant or susceptible than the resistant and susceptible parents, was observed in the RIL population for both the field and greenhouse experiments. Transgressive segregation has been reported in several studies (Jinks and Pooni 1976;Pilet-Nayel et al. 2002;Feng et al. 2011;Li et al. 2012;Coyne et al. 2015). The factors responsible for transgressive segregation of the progeny remain unclear (Kuczyn′ska et al. 2007). However, Nakedde et al. (2016) suggested that resistance genes in the parents residing on different linkage groups could account for the higher levels of tolerance exhibited by some of the RILs. The RILs exhibiting greater tolerance to ARR than the parents could be used in genetic crosses to stack the resistance genes, and the potential of the developed lines for various breeding programs should be exploited. Furthermore, sequence comparison between the tolerant parent and the RILs that showed greater tolerance could enable the identification of tightly linked markers for marker assisted breeding.
Although the G × E interaction for all traits were significant, moderate to high correlation coefficients were observed for all the traits between the field and greenhouse, as well as among individual experiments in the field and greenhouse and in the BLUPs data. High heritability of traits in each single field or greenhouse experiment, as well as in the BLUPs data, also confirmed the significant genetic effects on ARR tolerance in RILs. We observed a high correlation for disease severity and other traits in the field trials and greenhouse experiments (Table S1). Previous studies in chickpea (Cicer arietinum) (Johansen et al., 1994) and snap bean (Phaseolus vulgaris) (Navarro et al. 2008) indicated that early vigor had a beneficial effect on shoot biomass production. The epistatic (QTL × QTL) analysis showed a significant interaction of genomic regions linked to root rot severity, foliar weight and vigor (Table 4). The observation that about 81.5% (22 out of 27) of the QTL × QTL interactions were associated with root rot severity, foliar weight and vigor suggests that the same genomic regions control these traits. In the case of height, only 7.4% (2 out of 27) of the QTL × QTL interactions were associated with additive-effect markers for plant height. This suggested that plant height was a poor measure for ARR severity in pea. Thus, the results of our study are consistent with the findings of Conner et al. (2013), who suggested that ARR affected foliar weight and plant vigor but not plant height.
In the QTL mapping, genomic regions corresponding to 27 QTL for root rot severity, two disease-related criteria (foliar weight and vigor) and one agronomic trait (plant height) were identified. The largest (R 2 of 28.3-47.7%) major-effect QTL, Ae MRCD1 Ps-3.1, for resistance to A. euteiches identified in this study was found on LG III 1 3 (chrom 5), while the second largest (R 2 15.4-18.6%) majoreffect QTL, Ae MRCD1 Ps-4.1, was located on LG IV (chrom 4) ( Table 3). The QTL on LG IV (chrom 4) was the most stable, since it was detected in four of the five experiments, while the QTL on LG III (chrom 5) was detected in only one experiment (BLUPs data not counted). The largest majoreffect QTL detected by Pilet-Nayel et al. (2002, APh1 (R 2 of up to 47%), was also located on LG IV (chrom 4). Weeden (2000) also reported that a gene influencing tolerance to common root rot in pea was located on LG IV (chrom 4). In contrast, the third largest major-effect QTL (AePs4.1) detected by Hamon et al. (2011) was located on LG IV (chrom 4). The second major-effect QTL detected by Hamon et al. (2011) (Ae-Ps3.1; R 2 of up to 29.9%) was located on LG III (chrom 5). Hamon et al. (2011) detected Ae-Ps3.1 consistently in multiple experiments, while in this study Ae MRCD1 Ps-3.1 was identified in only one of the greenhouse experiments.
The QTL for resistance to A. euteiches isolate Ae-MRDC1 located on LG II (chrom 6), Ae MRCD1 Ps-2.1, was found to be a moderate-effect QTL with R 2 = 14.1 − 14.5%, while those on LG VII (chrom 7), Ae MRCD1 Ps7.1 and Ae MRCD1 Ps-7.2, were found to be minor-effect and moderate-effect QTL, respectively, with R 2 values ranging from 7.7 to 17.2%. The QTL reported on LG II (chrom 6) and LG VII (chrom 7) by Hamon et al. (2011) were found to be a combination of minor-effect and major-effect QTL. Ae-Ps2.2 was a majoreffect QTL (R 2 = 26.9%), while Ae-Ps2.1 and Ae-Ps2.3 were minor-effect QTL, accounting for up to 15.4% of the phenotypic variation. Similarly, two QTL Ae-Ps7.6a and Ae-Ps7.6b detected on LG VII (chrom 7) by Hamon et al. (2011) were moderate-effect (R 2 = 14.4%) and major-effect QTL (R 2 = 42.2%), respectively. Coincidentally, Pilet-Nayel et al. (2002 also detected two QTL (Aph6 and APh13) on LG VII (chrom 7) for disease-related criteria, namely aboveground index and root weight loss. Similar to the findings of Pilet-Nayel et al. (2005), the most stable and consistently detected QTL with the largest R 2 for aboveground disease indices (foliar weight and vigor) were located on LG IV (chrom 4), with R 2 values of up to 29.3% and 28.7%, respectively. Therefore, the major-effect and moderate-effect QTL detected in this study and those reported by Weeden et al (2000), Pilet-Nayel et al. (2002; and Hamon et al. (2011), appear to be on similar chromosomes, despite the different pedigrees of the parents.
The similarity in the number of major-and moderate-QTL detected in this and earlier studies may reflect the fact that commercial pea cultivars have been developed from a very limited pool of partially resistant pea germplasm (Lockwood and Ballard 1960;Shehata et al. 1983;Gritton 1990Gritton , 1995Kraft 1992;Davis et al. 1995;Wicker et al. , 2003Roux-Duparque et al. 2004;Pilet-Nayel et al. 2007). One of the progenitors of the tolerant parent '00-2067' used in this study was the plant introduction (PI) line 257593. PI 257593 was reported to be highly resistant to root rot caused by Fusarium and Pythium species (Kraft 1974). Only a handful of workers have focused on the development of pea cultivars partially resistant to root rot-causing pathogens including A. euteiches (Lockwood 1960;Kraft 1974Kraft , 1984Kraft , 2001Kraft and Burke 1974;Giles 1976, 1978;Gritton 1990;Davis et al. 1995;Gritton 1995;Kraft and Coffman 2000a, b;Roux-Duparque et al. 2004;Pilet-Nayel et al. 2002. Many breeding programs around the world have utilized the few partially resistant progenitor pea germplasm from North America. Some of the partially resistant germplasm might have also been crossed with each other to stack the resistance genes. Hence, it is likely that the partially resistant parent '00-2067' used in this study may share some common genetic basis with pea germplasm used in the previously reported studies. Marker-assisted selection requires the development of molecular markers either from the gene controlling the trait under study or from genomic regions flanking the gene. However, unlike qualitative traits involving dominant genes, MAS has not always been successful for quantitative traits involving polygenic genes (Xu and Crouch 2008;Hospital 2009). An obvious challenge is the stacking of the many genes controlling a complex trait into a single germplasm. In addition, many QTL are unstable in different environments, even if they have large effects, and negative epistatic interactions may reduce the efficiency of MAS (Hospital 2009). In this study, the peak genomic regions corresponding to Ae MRCD1  and Ae MRCD1 Ps-4.2, Fwt-Ps4.2 and Vig-Ps4.2 were within ≈ 3.0-5.0 cM of each other. The co-localization of root rot severity, foliar weight and vigor suggest that this region is very important for the resistance of pea to ARR. Therefore, targeting the entire 38.0-58.0 cM region on LG IV (chrom 4) may be an important breeding objective.
Based on our linkage map, the SNP marke r s P s C a m 0 3 7 5 4 9 _ 2 2 6 2 8 _ 1 6 4 2 a n d PsCam026054_14999_2864 bordered this region. The region contained 80 additional SNP markers (i.e., 4 markers/ cM). Markers in this region belonged to linkage disequilibrium (LD) block IV.8 in the study of Desgroux et al. (2016). Genes in LD block IV.8 included an LRR serine threonine protein kinase and vacuolar amino acid transporter, which are involved in the plant defense response, the FYVE zinc finger domain involved in signal transduction, AP2-like ethylene-responsive and bHLH123 transcription factors, and ATPase involved in biochemical and other cellular processes (Desgroux et al. 2016). Similarly, the genomic region corresponding to Ae MRCD1 Ps-3.1 and Vig-Ps3.1 on LG III (chrom 5), associated with partial resistance to ARR and vigor, respectively, is important for breeding. The peak region, however, was vast (5.0-33.0 cM) and contained fewer markers. As such, more markers need to be developed to fine map this genomic region as well as to identify candidate ARR resistance genes for future cloning and gene functional analysis. One approach could be to sequence the parents and each individual RIL, although genotyping costs may be prohibitively high for such a large number of genotypes. Alternatively, resequencing of important genomic regions, such the identified regions on LG IV (chrom 4) and LG III (chrom 5), would be more cost-effective. Recently, bulk segregant RNA-sequencing (BSR-seq) technology has emerged as a novel tool for the study of disease resistance and other traits in important crops (Liu et al. 2012;Yu et al. 2016;Hu et al. 2019;Wu et al. 2018a, b). BSR-seq in pea could help to identify the genes involved in various biological processes, regulation and development, as well as those involved in defense.
Previous studies evaluated partial resistance mainly by using disease severity (Pilet-Nayel et al. 2002Hamon et al. 2011Hamon et al. , 2013Lavaud et al. 2015;Desgroux et al. 2016) and parameters related to underground losses in pea (Pilet-Nayel et al. 2002Hamon et al. 2011Hamon et al. , 2013Desgroux et al. 2016). While the trait names and scoring scales were different, the vigor in our study was similar to above ground index (AGI) (Pilet-Nayel et al. 2002 or aerial decline index (ADI) (Hamon et al. 2011(Hamon et al. , 2013Desgroux et al. 2016). Pilet-Nayel et al. (2002 also used the trait percentage of dried weight losses (DWL), while we applied dry foliar weight to indicate the effect of ARR on aboveground biomass. Thus, the parameters measured in this study were consistent with those used in the aforementioned studies.
In conclusion, linkage analysis and QTL mapping using high-density SNP markers and SSR anchor markers and an F 8 RIL population enabled us to identify a 20.0 cM chromosomal region on LG IV (chrom 4) and a 28.0 cM genomic region on LG III (chrom 5) as being largely responsible for partial resistance to A. euteiches isolate Ae-MDCR1. Extensive validation of the identified markers is needed to determine their utility, particularly given the challenges associated with MAS of quantitative traits.
Acknowledgements The Canadian Agricultural Partnership (CAP) provided financial support for this research under project # 1000210132. The University of Alberta and Alberta Agriculture and Forestry (Crop Diversification Centre North) provided in-kind support, including access to greenhouse facilities and molecular laboratory equipment. Waldo C. Penner and Dennis B. Stoesz provided technical help, while Dr. Kenneth B. McRae helped in field designs and provided statistical advice on the analysis of the field research results.
Author contribution LFW contributed to isolation of Aphanomyces euteiches isolate Ae-MRDC1, inoculum preparation, greenhouse screening of RIL population and parents for resistance to Aphanomyces root rot, disease rating, measurement of foliar weight, vigor and plant height, phenotypic data analysis, DNA extraction, PCR, genotyping with SSR markers and writing of the manuscript. RFA contributed to supervision of the molecular marker work, molecular data analysis, linkage map construction, QTL mapping and writing of the manuscript. SFH was the principal investigator and contributed to grant application, supervision and provision of technical support to the graduate student for RIL population screening in the greenhouse and revision of the manuscript. RLC contributed to development of the RIL population, field evaluation of RIL, disease evaluation, measurement of vigor and foliar weight and revision of the manuscript. KFC contributed to grant application, supervision and provision of technical support to the graduate student for RIL population screening in the greenhouse. DLM contributed to grant application and the provision of technical support to the graduate student. SES was the principal investigator and contributed to grant application, supervision and provision of technical support to the graduate student and revision of the manuscript.
Funding Institutions supporting this project are indicated in the acknowledgements section.
Availability of data and material Heritability of each trait, marker information, linkage information and QTL profiles are available in the main manuscript or as supplementary data.

Declarations
Conflict of interest On behalf of all authors, the corresponding author states that there is no conflict of interest.
Ethical approval On behalf of all authors, the corresponding author states that all experiments in this study complied with the ethical standards in Canada.