A High-Density Integrated DArTseq SNP-Based Genetic Map of Pisum fulvum and Identification of QTLs Controlling Rust Resistance

Pisum fulvum, a wild relative of pea is an important source of allelic diversity to improve the genetic resistance of cultivated species against fungal diseases of economic importance like the pea rust caused by Uromyces pisi. To unravel the genetic control underlying resistance to this fungal disease, a recombinant inbred line (RIL) population was generated from a cross between two P. fulvum accessions, IFPI3260 and IFPI3251, and genotyped using Diversity Arrays Technology. A total of 9,569 high-quality DArT-Seq and 8,514 SNPs markers were generated. Finally, a total of 12,058 markers were assembled into seven linkage groups, equivalent to the number of haploid chromosomes of P. fulvum and P. sativum. The newly constructed integrated genetic linkage map of P. fulvum covered an accumulated distance of 1,877.45 cM, an average density of 1.19 markers cM−1 and an average distance between adjacent markers of 1.85 cM. The composite interval mapping revealed three QTLs distributed over two linkage groups that were associated with the percentage of rust disease severity (DS%). QTLs UpDSII and UpDSIV were located in the LGs II and IV respectively and were consistently identified both in adult plants over 3 years at the field (Córdoba, Spain) and in seedling plants under controlled conditions. Whenever they were detected, their contribution to the total phenotypic variance varied between 19.8 and 29.2. A third QTL (UpDSIV.2) was also located in the LGIVand was environmentally specific as was only detected for DS % in seedlings under controlled conditions. It accounted more than 14% of the phenotypic variation studied. Taking together the data obtained in the study, it could be concluded that the expression of resistance to fungal diseases in P. fulvum originates from the resistant parent IFPI3260.


INTRODUCTION
The cultivated pea (Pisum sativum L. subsp. sativum) is one of oldest domesticated crops. Since its domestication about 10,000 years ago it has been improved for important agronomic traits, being today the cool season grain legume most cultivated in Europe and the second in the world (Rubiales et al., 2015). However, pea yield is still relatively unstable due to its limited adaptability to a wide range of environmental conditions and its susceptibility to diseases and pests. Among the biotic stresses, the rusts have been acquiring more and more economic importance since the mideighties. Uromyces pisi (Pers.) Wint. is the main causal agent of pea rust in temperate regions (Barilli et al., 2009a,b) although in tropical and subtropical regions as India or China, it has been reported to be caused by the fungus Uromyces viciae-fabae (Pers.) J. Schröt [syn. U. fabae (Pers.) de Bary] (Singh and Sokhi, 1980;Kushwaha et al., 2006). In years of epidemic the rust causes the drying of leaves and severely affects the development of pods, with consequent yield losses that can reach up to 50% (Barilli et al., 2009a).
Unfortunately, there is little resistance available in the cultivated P. sativum, being insufficient to achieve an effective control. Only some levels of incomplete resistance have been identified after an exhaustive screening of U. pisi resistance in a large germplasm collection identifying the highest levels of resistance in P. fulvum Sibth. & Sm. (Barilli et al., 2009a,b,c). Not without difficulty, crosses between P. sativum and P. fulvum have been possible and, in fact, resistance genes to Bruchus pisorum and Erysiphe pisi have already been transferred from P. fulvum to P. sativum by natural hybridization (Fondevilla et al., 2007;Byrne et al., 2008). Similarly, resistance to other biotic stresses such as Didymella pinodes or Orobanche crenata has been identified in P. fulvum, and crosses have also been carried out to transfer the resistance genes to the cultivated pea (Rubiales et al., 2005(Rubiales et al., , 2009Fondevilla et al., 2008Fondevilla et al., , 2017. The use of modern breeding tools will facilitate the efficient exploitation of all the potential offered by P. fulvum for P. sativum breeding. DNA-based genetic markers provide powerful tools for the identification and localization of genes that control traits of agronomic importance, whose subsequent introgression into commercial varieties facilitates the advancement of plant breeding programs. Despite the important source of resistance provided by P. fulvum so far there has been relatively little research on the development and application of molecular genetic tools in wild peas compared to the best studied P. sativum in which several types of markers have already been related to several disease resistance regions (Pavan et al., 2013;Carrillo et al., 2014;Coyne et al., 2015;Sudheesh et al., 2015;Tayeh et al., 2015;Boutet et al., 2016).
The association of molecular markers with U. pisi resistance has been previously studied in the F 2 population derived from the intraspecific cross between the P. fulvum accessions IFPI3260 × IFPI3251 (Barilli et al., 2010). A previous genetic map of 1283.3 cM that included 146 markers (144 RAPDs and 2 STSs markers) distributed in 9 linkage groups (LGs) was developed. It was observed that rust resistance was governed by a major QTL that explained 63% of the total phenotypic variation located in the LG 3 and flanked by the RAPD markers OPY11 1316 and OPV17 1078 . To our knowledge, there is no other report on QTLs for U. pisi resistance in pea. Vijayalakshmi et al. (2005) and Rai et al. (2011) studied in a BC 1 F 2 population the resistance to U. fabae Pers. de-Bary, the pathogen responsible for rust disease in pea in sub-tropical regions, reporting that it was governed by a major and a minor QTL identified in the LG 7 as Qruf (22.4-58.8% of phenotypic variation; flanked by SSR markers AA446 and AA505) and Qruf1 (11.2-12.4% of phenotypic variation; flanked by SSR markers AD146 and AA416).
As mentioned above, several pea linkage maps based on different type of markers have been constructed. Most of these are gel-based markers and have limited ability to rapidly analyse a large number of marker loci. Although some of these limitations can be overcome by using specialized hardware such as highthroughput capillary electrophoresis equipment that can improve the ability of allelic discrimination, reproducibility, and speed, there are still limitations related to the sequential nature, the requirement to have DNA sequence information to expand the currently available marker toolkits and the high costs of use these marker technologies. Mainly this last factor has made the use of molecular markers directly in wild species impractical. In this scenario, the Diversity Arrays Technology (DArT) in combination with next-generation sequencing platforms (Kilian et al., 2012;Raman et al., 2014) known as DArTseq TM , provides a good choice as a high throughput marker genotyping platform that can develop a relatively large number of polymorphic markers to build dense genetic maps with low-cost investments. The thorough coverage of the genome and the high-density genetic maps based on DArTseq TM technology increase the power of QTL detection (Thudi et al., 2011). Additional advantages of DArTseq TM technology are its suitability for polyploid species as well as the possibility of developing rapidly for virtually any genome. As a result, DArTseq-derived markers are currently used in more than 400 species (http://www.diversityarrays. com/), being very popular among crops with the non-sequenced genomes. The DArTseq markers include genome-wide profiling of a large number of SNPs and detection of insertion/deletion polymorphisms which is easily expandable for genetic scope (Kilian et al., 2012). So far, the DArT markers have been used in legumes for the mapping of pigeon peas (Cajanus cajan L.) (Yang et al., 2011), chickpeas (Cicer arietinum L.) (Thudi et al., 2011), lupins (Lupinus albus L.) (Vipin et al., 2013), common beans (Phaseolus vulgaris L.) (Valdisser et al., 2017), and soybeans (Glycine max L.) (Vu et al., 2015).
Therefore, the goals of this work have been the development of the first integrated high-density genetic linkage map of P. fulvum and the identification of QTLs of resistance to pea rust.

Plant Material
A mapping population of 84 F 7 RILs derived from the cross IFPI3260 (resistant) × IFPI3251 (susceptible) was used to construct an integrated SilicoDArT + SNP + SSR + STS based linkage map and was screened for rust (U. pisi) resistance.

Phenotyping Analysis for Rust Resistance under Field and Controlled Conditions
Seeds of the pea susceptible check cv. Messire and 84 F 7 RILs from the cross P. fulvum IFPI3260 (resistant) × IFPI3251 (susceptible) segregating for rust resistance and their parental lines were sown in fields located at Córdoba (Spain) during the 2013/2014, 2014/2015 and 2015/2016 growing seasons. Each accession was represented by a single row of 1 m in which 15 plants were sown and which was separated from the adjacent row by 0.7 m. The assays for each season included three replications of each accession arranged in a completely randomized design. Plots were artificially inoculated at the end of March to ensure high and uniform levels of rust infection. Prior to field inoculation, stock spores of the local U. pisi isolate Up-Co01 conserved at −80 • C were heat shocked at 40 • C for 5 min and multiplied on susceptible pea cv. Messire. These seedlings were sprayinoculated with an aqueous urediniospores suspension (±1.0 × 10 5 urediniospores ml −1 ) to which Tween-20 (0.03%, v: v) was added as a wetting agent and incubated for 24 h at 20 • C in complete darkness and 100% relative humidity. Later, they were transferred to a growth chamber at 20 • C with a photoperiod of 14 h of light and 10 h of darkness and a light intensity of 148 µmol m −2 s −1 . After 2 weeks, the fresh urediospores were collected using a vacuum spore collection device and used for field inoculation. The plants were inoculated after sunset to benefit from the darkness and high relative humidity of the night. At maturity, disease severity (DS) and infection type (IT) were assessed. DS was visually estimated as the percentage of canopy covered by rust pustules. IT was assessed using the 0-4 scale of Stackman et al. (1962), where IT 0 = no symptoms, IT; = necrotic flecks, IT 1 = minute pustules barely sporulating; IT 2 = necrotic halo surrounding small pustules, IT 3 = chlorotic halo and IT 4 = well-formed pustules with no associated chlorosis or necrosis.
In addition, the response of the plant material was studied in seedlings under controlled conditions. Two seeds were sown in each pot (35 × 35 cm) with a 1:1 mixture of sand and peat. Three replicates were performed with 8 plants (4 pots) per accession and replicate in a completely randomized design. Seedlings were inoculated when the third leaf was completely expanded (± 12 days after germination). The inoculation was carried out by dusting the plants with the urediospores of U. pisi isolate Up-Co01 (2 mg spores pot −1 ) diluted in pure talc (1:10, v:v) using a spore settling tower. The plants were incubated for 24 h at 20 • C in complete darkness and at 100% relative humidity. Then they were transferred to a growth chamber at 20 • C with a photoperiod of 14 h of light and 10 h of darkness and 148 µmol m −2 s −1 of irradiance at plant canopy. DS and IT values were recorded 10 days after the inoculation as above mentioned. DS values measured in both field and controlled conditions were normalized by the arcsine transformation [y = arcsine ( √ DS)].
Pearson's linear correlations between rust resistance parameters were performed using Statistix (version 8.0; Analytical Software, Tallahassee, USA).

DNA Extraction and Quantification
The eighty-four P. fulvum F 7 RILs of the mapping population and the parents of the cross were grown under controlled conditions at IAS-CSIC. Around 1 g of the young leaf tissue from the 3rd to 4th node of each seedling was excised, immediately frozen in liquid nitrogen and stored at −80 • C. Genomic DNA was isolated from the frozen leaves using a modified cetyltrimethylammonium bromide (CTAB)/chloroform/isoamylalcohol method (Doyle and Doyle, 1987). DNA quantification was performed by agarose gel electrophoresis (0.8 %), and it was adjusted to 50 ng/µl for DArT and SNP genotyping and to 5 ng/µl for SSR genotyping.

Genotyping of Individual DNA Samples Using DArTseq Technology
A high-throughput genotyping method using the DArT-Seq TM technology at Diversity Arrays Technology Pty Ltd (Canberra, Australia) was implemented to genotype the F 7 RIL population of P. fulvum. Essentially, DArT-Seq TM technology relies on a complexity reduction method to enrich genomic representations with single copy sequences and subsequently perform nextgeneration sequencing using HiSeq2000 (Illumina, USA). DArT-Seq detects both SNPs and presence-absence sequence variants, collectively referred to as DArT-Seq markers (Raman et al., 2014). DArT-Seq was optimized for P. fulvum by selecting the most appropriate complexity reduction method (PstI-MseI restriction enzymes). DNA samples were processed in digestion/ligation reactions as described by Kilian et al. (2012), but replacing a single PstI-compatible adapter with two different adapters corresponding to two different restriction enzymes (RE) overhangs. The PstI-compatible adapter was designed to include Illumina flowcell attachment sequence, sequencing primer sequence, and staggered, varying length barcode region. The reverse adapter contained the flowcell attachment region and MseI-compatible overhang sequence. Only "mixed fragments" (PstI-MseI) were effectively amplified in 30 rounds of PCR using the following reaction conditions: 1 min at 94 • C for initial denaturation; 30 cycles each consisting of 20 s at 94 • C for denaturation, 30 s at 58 • C for annealing and 45 s at 72 • C for extension; and finally a 7 min extension step at 72 • C. After PCR, equimolar amounts of amplification products from each sample of the 96-well microtiter plate were bulked and applied to c-Bot (Illumina) bridge PCR followed by sequencing on Illumina Hiseq2000. The sequencing (single read) was run for 77 cycles. Sequences generated from each lane were processed using proprietary DArT analytical pipelines. In the primary pipeline, the FASTQ files were first processed to filter poorquality sequences, applying more stringent selection criteria to the barcode region compared to the rest of the sequence. Thus, the assignments of the sequences to specific samples carried in the "barcode split" step were more consistent. Approximately 2,500,000 (±7%) sequences per barcode/sample were used in marker calling. Finally, identical sequences were collapsed into "fastqcall files." These files were used in the secondary pipeline for DArT P/L's proprietary SNP and SilicoDArT (Presence/Absence Markers in genomic representations) (present = 1 vs. absent = 0) calling algorithms (DArTsoft14). The analytical pipeline processed the sequence data. Parameters for the DArT marker assaying pipeline for quality control (Kilian et al., 2012) such as: (i) the reproducibility of 100%, (ii) the overall call rate (percentage of valid scores in all possible scores for a marker) over 95%, (iii) the polymorphic information content (PIC) between 0.3 and 0.5 and (iv) the Qvalue (the logarithm of the minimum false discovery rate at which the test may be called significant) above 2.5 were used for selecting high-quality SilicoDArT and derived SNPs markers for genetic mapping.

Genotyping with Simple Sequence Repeats (SSRs), Sequence Tagged Site (STS) and Single-Nucleotide Polymorphism (SNP) Markers
A set of 46 genic and genomic SSRs (15 fluorescently labeled SSRs and 31 M13-labeled SSRs) previously described by Loridon et al. (2005) was screened using the two parental lines and the 84 F 7 RILs. A multiplex PCR (markers labeled with different dyes) was carried out to amplify SSRs following optimal PCR conditions previously published by Loridon et al. (2005). Then 2 µl of the PCR product were taken from each marker of the multiplex set and pooled together for simultaneous detection of the amplified alleles. Seven µl of formamide and 0.2 µl of fragment-size standard GeneScanTM 500 LIZ were added to the pooled PCR product and run on an ABI 3730 DNA genetic analyzer (Applied Biosystems). The data were collected automatically by the detection of the different fluorescences and analyzed using GeneMapper v4.0 software (Applied Biosystems).
In addition 10 STS markers reported by Gilpin et al. (1997) were surveyed for polymorphism using the protocol described by Barilli et al. (2010). When no polymorphism was detected, PCR products from both parents were digested with a range of restriction endonucleases which recognized 4-and 5-base sequences (BioLabs inc , Ipswich, MA, USA). A 0.2 µl aliquot of restriction enzyme, 12 µl of sterile water, and 2.5 µl of the buffer required for each enzyme were added to 10 µl of the PCR reaction, and the digestion was incubated overnight at 37 • C. View of PCR products and gel images analyses were performed as described by Barilli et al. (2010).

Linkage Mapping and QTL Mapping
The scores of all polymorphic DArTseq markers, SSR, STS, and SNP markers were converted into genotype codes ("A, " "B") according to the scores of the parents. Linkage groups were obtained using the software JoinMap version 4.1 (Van Ooijen, 2006). The maximum likelihood mapping algorithm, which was optimized for constructing dense genetic maps using this software (Jansen et al., 2001), was first used for grouping all of the polymorphic markers. Then, the method of regression mapping (Haley and Knott, 1992) was used for map construction with approximately 1,000 markers with appropriate genetic distance and the marker position and the order of markers for three rounds to merge the tightly adjacent markers into bins. The markers in adjacent loci with genetic distance below 0.2 cM were classified into a bin during the first two rounds of mapping. Moreover, one marker with sequence information and with the least missing genotype from each bin was chosen as a "bin representative" for the next round of genetic mapping. For the last (the third) round of mapping, the makers in adjacent loci pairs with genetic distances below 0.1 cM were classified into a bin to avoid incorrect classification when the markers were decreased in the map. The Kosambi mapping function (Kosambi, 1943) was used to convert recombination frequencies into map distances, and only "Map 1" was used for further analysis. The LG maps of each chromosome were drawn and aligned using MapChart v2.3 (Voorrips, 2002). The threshold for the goodness of fit was set to ≤ 5.0, with a recombination frequency of <0.4. Linkage groups were separated using a LOD score ≥ 3.0. Markers with a mean Chi-Squared value of recombination frequency > 4.0 were discarded. The DArT markers were named with the numbers corresponding to unique clone ID following Kilian et al. (2012).
Several previously described anchor markers were used to find the correspondence between P. fulvum and P. sativum linkage groups and assign them to pea chromosomes. For the same purpose, the sequences from DArT-seq-derived markers were compared with Medicago truncatula genomic backbone by using Phytozome v.12 (https://phytozome.jgi.doe.gov/pz/portal.html) to perform a synteny analysis using three parameters recently defined by Salse et al. (2009). These parameters increase the stringency and significance of BLAST sequence alignment by parsing BLASTX results and rebuilding HSPs (High Scoring Pairs) or pairwise sequence alignments to identify accurate paralogous and orthologous relationships. This analysis allowed searching for sequence similarity-based homology between legume species providing an alternative approach to finding correspondence between linkage groups. QTL analysis for rust resistance was conducted using composite interval mapping (CIM) and multiple interval mapping (MIM) in MapQTL 6.0 package (Van Ooijen, 2011). Markers to be used as cofactors for CIM were selected by forward-backward stepwise regression. Significance thresholds of log of odds (LOD) corresponding to a genome-wide confidence level of P < 0.05 were determined for each trait using the permutation test of MapQTL 6.0 with 1000 iterations, according to Barilli et al. (2010). Skewness and Kurtosis coefficients were calculated following Lynch and Walsh (1997). The coefficient of determination (R 2 ) for the marker most tightly linked to a QTL was used to estimate the proportion of the total phenotypic variation explained by the QTL.

Rust Resistance
Pisum fulvum accessions showed compatible interaction (high IT) against U. pisi although they differed in DS, confirming previous findings (Barilli et al., 2009a,b). Analysis of variance of each trial, either under field or controlled conditions revealed highly significant genotypic effects for DS resistance criteria among the RIL families (P < 0.05). The parental accessions showed contrasting DS responses, with IFPI3260 being highly resistant (DS < 5%) and IFPI3251 highly susceptible (DS > 30% in field and DS = 63% in controlled conditions) (Figure 1). Distribution of residuals after analysis of variance was normal in each year and tested conditions according to Shapiro-Wilk (P > 0.01). Variances of genotypes and replicates were homogeneous according to Bartlett's test (P > 0.05). The coefficient of Skewness in field conditions was of 0.86, 0.87, and 0.83 for seasons 2013/14, 2014/15, and 2015/16, respectively, which indicated that the population distribution tended to the resistance as the parent IFPI3260. Distributions of the estimated adjusted means are represented in Figure 1 for DS for each year and condition evaluated. DS values ranked between 0.75 and 45% under field conditions and 5-63.3% under controlled conditions. They did not differ from the normal distribution, confirming the quantitative inheritance of the partial resistance. In addition, transgressive segregants with increased resistance and susceptibility compared with the parentals were observed for DS resistance criteria over years and conditions. Pearson's linear correlation coefficients (Table 1) were highly significant between DS values evaluated within years in the field, as well as under controlled conditions.

Genotyping and Linkage Mapping
A total of 9,569 high-quality SilicoDArT and 8,514 SNPs markers were identified. Of these, a set of 8,538 SilicoDArT markers (89.2%) and 3,483 SNPs markers (40.9%) were selected for mapping after quality filtering. Thirty-five of the set of 46 SSRs markers previously described for P. sativum by Loridon et al. (2005) amplified in P. fulvum, which corresponded to a potential transferability rate of 76% for this marker type. However, only 12 markers were polymorphic. This meant that only 26.1% of the SSR markers tested showed utility for mapping in our population. The remaining 21 yielded non-specific amplification and were excluded from further analysis. The mapping dataset was complemented with 20 of the 26 SNP markers (76.9%) developed using the BeadXpress Primer Design (Illumina, San Diego, CA, USA) (Deulvot et al., 2010), as well as with 5 of the 10 STS markers (50%) reported by Gilpin et al. (1997), which provided reliably scorable polymorphism. As a result, 37 non-DArT-based markers were added to a set of 12,021 DArTseq   derived markers that scored co-dominantely (Tables 2, 3). All the non-DArT-based markers mapped in the expected LGs according to previous publications (Loridon et al., 2005;Deulvot et al., 2010;Bordat et al., 2011;Carrillo et al., 2014;Tayeh et al., 2015). Markers were distributed across 7 LGs using LOD thresholds ranging from 3 to 10 and a recombination frequency (r) threshold < 0.4 (JoinMap vs. 4) (Supplementary File 1). Each assigned group included at least two markers common to other published P. sativum genetic maps. The newly constructed integrated genetic linkage map of P. fulvum covered a total length of 1877.45 cM, with an average density of 1.19 markers cM −1 and an average adjacent-marker gap distance of 1.85 cM ( Table 2). The total number of mapped loci per linkage group ranged from 1347 on LGVI to 2068 on the LGII, and the average was 1,723 loci LG −1 . The longest individual linkage group map was for the LGII (349.48 cM), the shortest was for the LGIV (219.54 cM) (Figure 2), and the average LG length was 268.21 cM. The density of markers in the individual linkage groups ranged from 1.76 markers cM −1 in the LGVII to 2.09 markers cM −1 in the LGVI. Map distances between two consecutive markers varied from 0 to 13.04 cM, while the gap average between markers varied from 1.66 cM in the LGVII and 2.05 cM in the LGVI ( Table 2). Thirty-seven previously described markers (Table 3), as well as 98 sequences from P. fulvum DArTseq-derived markers that were BLASted with M. truncatula genome (Supplementary File 2), allowed to define the correspondence between P. fulvum and P. sativum linkage groups and the pea chromosome assignment, as follows: 2 SSR markers (AA317 and AA160) as well as 20 DArTseq-derived markers linked the P. fulvum LGI to the P. sativum LG7; 5 previously reported markers, as well as 14 DArTseq-derived markers linked LGII to LG5; 16 previously reported markers, as well as 12 DArTseq-derived markers linked LGIII to LG3; 4 markers (P482-HhaI, AA102, peptrans_SNP1 and peptrans_SNP2) and 5 DArTseq-derived markers related LGIV to LG2; 4 previously reported markers (AA335, ARBA3199_567, PsDHN1_124, and PsDHN1_320) and 10 DArTseq-derived markers related the LGV to the P. sativum LG6; 3 SNP markers (agps2_SNP3, ARBA8726_644, and subt_SNP2) and 18 DArTseq-derived markers linked the LGVI to the P. sativum LG1; finally, 3 previously mapped markers (P393-HhaI, GTPaseact_SNP1 and gpt2_SNP1) and 10 DArTseq-derived markers showed that the LGVII corresponds with the LG4. The LGVI displayed an inverted markers order compared to Bordat et al. (2011) and Carrillo et al. (2014). Pea chromosome assignment following Ellis and Poyser (2002) was reported in Table 3.

QTL Mapping
The scoring of disease severity under controlled and field conditions allowed the identification of genomic regions involved in pea rust resistance. Quantitative trait loci analysis using CIM and MIM methods revealed 2 genomic regions associated with the resistance of adult plants to U. pisi under field conditions located in the LGs II and IV. Results were consistent along all the growing seasons considered. The QTL UpDSII showed a LOD score range of 3.28-4.51 and explained 20.7-29.2% of rust severity variation at the adult plant stage (Table 4, Figure 3A). UpDSII was localized between 163.7 and 164.2 cM from the beginning of the LGII, between the derived DArT marker 3567800 and the previously known SSR marker AD280. The distance to the left and to the right flanking markers ranged between 0.9-1.4 and 1.2-1.7 cM, respectively. The second consistent QTL, UpDSIV, showed a LOD score ranging from 3.00 to 3.47 and explained 19.8-22.1% of the variation of the rust severity at the adult plant stage (Table 4, Figure 3B). UpDSIV was localized between 25.4 and 26.7 cM from the beginning of the LGIV, between the derived DArT markers 3563695 and 3569323 (for DSfie14 and DSfie15), and 3560101_51:A>T and 3536169 (for DSfie16) ( Table 4). The distance to the left and to the right flanking markers ranged between 1.5-4.7 and 0.9-1.5 cM, respectively. Both QTLs explained together 43.8, 40.5, and 51.3% of pea rust severity variation in DSfie14, DSfie15, and DSfie16, respectively.
These results were also confirmed under controlled conditions after performing the QTL analysis for resistance to U. pisi in wild pea seedlings. In fact, the consistent QTL UpDSII (LOD = 4.42) located in the LGII was localized at 165.4 cM from the beginning of the LGII (peak QTL position coinciding with the SSR marker AD280), between the derived DArT markers 3534625 and 3539148. The distance to the left and the right flanking markers was of 1.3 and 1.5 cM, respectively ( Figure 3A). The second consistent QTL, UpDSIV, was also located at 25.4 cM from the beginning of the LGIV (LOD = 3.05) between the derived DArT markers 3563695 and 3569323, as previously found ( Figure 3B). The distance to the left and to the right flanking markers was of 4.2 and 1.4 cM, respectively and explained 19.9% of the phenotypic variation (Table 4). In addition, a third QTL associated with rust resistance in controlled conditions was found in the LGIV (named UpDSIV.2) ( Table 4). UpDSIV.2 showed a LOD score of 3.0 and explained 14.1% of the phenotypic variation. The distance to the left and to the right flanking markers was of 2.8-0.2 cM, respectively. All the three QTLs explained together 62.2% of pea rust severity variation in DSchamb (Table 4). Sequence information of all flanking markers is listed in Supplementary File 1.
The resistance-enhancing allele originates from the resistant parent IFPI3260 as shown by the negative value of the additive genetic effect (ranging from −13.11 to −2.93) ( Table 4).

DISCUSSION
Genetic linkage maps are highly valuable tools for the identification of candidate genes/QTLs that can be used in map-based cloning and marker-assisted breeding programs. Dense genetic maps based on sequence-derived markers, such as diversity arrays technology (DArTseq TM ) markers accelerate the process of fine mapping and map-based cloning of genes/QTLs. DArT-derived markers have been extensively used in recent years to build comprehensive genetic maps in legumes such as chickpea, pigeonpea, peanut and soybean (Vu et al., 2015;Ahmad et al., 2016). To the best of our knowledge, the present study is the first extensive published linkage map reported in wild pea (P. fulvum). The RIL mapping population was generated from a cross between two parental lines belonging to the P. fulvum specie showing excellent genetic diversity and providing a high frequency of polymorphisms of great utility for the construction of the linkage map. In fact, a total of 12,021 (66.5%) markers were found polymorphic in the panel of 84 RILs. In addition, the selected wild pea DArTseq markers showed a good average polymorphic information content (PIC) value of 0.42 (values ranging from 0.30 to 0.50), which indicated that those markers should be considered useful or informative (Kilian et al., 2012). Therefore it is very evident that DArT markers can be developed and typed quickly and economically compared to other marker technologies (Ahmad et al., 2016). The map of P. fulvum constructed in this work is composed by 7 LGs that cover a total length of 1,877.4 cM and an FIGURE 2 | Distribution of the DArT-seq-based and no-DArT-seq markers within each linkage group (LG) forming the map derived from the cross between P. fulvum IFPI3260 × IFPI3251. The x axis shows the LG and the y axis shows the genetic distance (cM).
average adjacent-marker distance of 1.8 cM. It includes 12,021 polymorphic DArT markers that together with 37 previously published markers of diverse nature (SNP, SSR and STS working as anchor markers) has been used as a bridge allowing us to combine our linkage map with recently published P. sativum consensus maps (Loridon et al., 2005;Bordat et al., 2011;Carrillo et al., 2014;Tayeh et al., 2015) and determine the orientation of the 7 LGs found. The comparison of these bridge markers of our IFPI3260 × IFPI3251 P. fulvum map with the previous P. sativum maps has revealed a general high consistency of marker order with only minor exceptions, as it was the change of the order of some marker within the same LG, or the case of the LGVI that showed an inverted order of anchoring markers with respect to those previously reported in P. sativum (Loridon et al., 2005;Carrillo et al., 2014). Such differences in the order of the markers among P. fulvum and P. sativum linkage maps were not unexpected, since genetic mapping only provides an indication of the relative positions and genetic distances of the markers among themselves (Loridon et al., 2005) and structural rearrangements of the chromosomes could be relatively common between P. fulvum and P. sativum (De Martino et al., 2000). Nonetheless, the inconsistency of the position on the map of these few SSRs could be explained by the presence of closely linked DArT loci (Kilian et al., 2012;Lan et al., 2017).
Over 85% of the 12,061 mapped markers showed markermarker linkage tendency (in groups of at least two markers), which collapsed into 994 bins. These markers were likely to represent mostly gene-rich regions, as DArT method of complexity reduction targets the hypomethylated regions of the genome. This is consistent with the observations found in the genetic maps of many other species such as chickpea and rapeseed (Thudi et al., 2011;Raman et al., 2014) as an example. Mapping of DArT array markers to the Eucalyptus reference genome using the unique sequence tag of each marker has suggested that PstI-based DArT markers are predominantly at the low-copy, gene-rich regions (Petroli et al., 2012). However, most of the observed clustering of marker could be explained by the limited resolution of our mapping population because its size was small (below 100 individuals).
The total genetic length of the presented map is 46.3% larger than the map (1,283.3 cM) previously constructed using the F 2 IFPI3260 × IFPI3251 RIL population in which 146 markers including 144 RAPD and 2 STS (Barilli et al., 2010) were mapped, but smaller than the maps of P. sativum published so far. In fact, the total lengths of some of them has even reached 2,416 (Timmerman-Vaughan et al., 2004) or 2,555 cM (Sudheesh et al., 2015). The map length obtained in the present study suggests that high-density and high-quality DArT genotyping data obtained after minimizing genotyping errors and missing data have contributed greatly to a better estimation of the distance among the markers.
The distribution of the markers was reasonably uniform along the map. In fact, the marker density of each individual linkage group ranged from 1 marker/1.764 cM (LGVII) to 1 marker/2.087cM (LGII), with a total average of 1 marker/1.89 cM with is a 9-fold increase compared to the previous published P. fulvum map (Barilli et al., 2010). It has been suggested that uniformly distributed loci every 10 cM over the entire genome is effective for MAS and QTL identification (Stuber et al., 1999). In this contest, the high level of genome coverage achieved in this map will be particularly useful for selecting markers for use in whole-genome breeding strategies and for saturating 4 | Position and effects of the quantitative trait loci (QTLs) for resistance to U. pisi in P. fulvum IFPI3260 × IFPI3251 RIL population based on the percentage of disease severity (DS%) scored in the field during three growing seasons (2013/14, 2014/15, 2015/16)   genomic regions of interest in other mapping populations. In addition, although around 12,000 markers were identified in this study, more markers are available for further work since more than 18,000 clones from P. fulvum were obtained. Thus, the P. fulvum DArT platform has been shown to be useful for the application in genome-wide screening for QTL discovery. It could also be expected to demonstrate its usefulness for recurrent parent background recovery in marker assisted backcrossing, for isolation of genes via map based cloning, comparative mapping and genome organization studies. The availability of better saturated molecular maps that can be achievabled using DArTSeq technology, as an example, will certainly provide breeders and geneticists with a much-desired tool for identifying genomic regions of interest, which in turn will increase the efficiency of marker-assisted breeding according to Ren et al. (2015). The present work is also the first report of a QTL analysis for phenotypic traits on a high-density linkage map in a F 7 derived RIL population of P. fulvum using the DArTseq TM technology. Pea rust caused by U. pisi is a major challenge to pea growers in temperate regions but complete resistance to this disease has not been identified so far (Barilli et al., 2009b,c). However, high incomplete resistance has been reported in wild Pisum accessions, especially in P. fulvum accession IFPI3260 that shows good levels of resistance at both seedling and adult plant stages and to different isolates of U. pisi (Barilli et al., 2012), being a suitable source of resistance to this disease. In the present study, no hypersensitive resistance reaction to U. pisi was observed. Neither the resistant parent nor the RILs were free of rust infection in the 3-year field experiments at Córdoba (Spain) and neither under controlled conditions. The resistance traits scored showed a continuous distribution in the RIL population, which indicated the quantitative nature of their inheritance. This is in agreement with the previous findings in which partial resistance was not associated with hypersensitive response against U. pisi. (Barilli et al., 2009a,b,c) or U. fabae (Singh and Sokhi, 1980;Vijayalakshmi et al., 2005;Chand et al., 2006;Ren et al., 2015) in Pisum spp. Most of the RIL families revealing a level of resistance skewed toward low disease severity rates which suggested that the combination of both parental lines could enhance the level of resistance provided by the resistant accession IFPI3260.
Two major and a minor QTLs distributed over two linkage groups were found to be associated with the DS%. Two consistent QTLs, UpDSII and UpDSIV, were located in the LGs II and IV, respectively. A minor QTL named UpDSIV.2 was also located in the LGIV. This was not in agreement with the previous QTL analysis performed on the F2 RIL population derived from the same cross, where a single genomic region associated with the resistance to U. pisi was identified and a QTL (Up1) explaining more than 60% of phenotypic variation under controlled conditions was described on LG 3 (Barilli et al., 2010). This discrepancy between both analyses could presumably be due to the advancement in the generation of the segregating populations, the number and nature of markers involved and/or software processing. The previous QTL Up1 was located in a genomic region poorly saturated in markers and the RAPD markers associated were not robust and transferable (Barilli et al., 2010).
The present study has demonstrated that the genetics of pea rust resistance is somewhat complex and controlled by several QTLs. The QTLs UpDSII and UpDSIV were consistently identified by DS% in both adult plants over 3 years at the field and in seedling plants under controlled conditions. Whenever they were detected, their contribution to the total phenotypic variance was substantially high (19.8-29.2). Previous studies in legumes have shown good correlations between rust disease assessments in seedlings under controlled conditions and adult plants in the field, suggesting the existence of genetic factors that control the effectiveness of the resistance at different developmental stages and environments (Singh et al., 2015;Rai et al., 2016). For U. fabae resistance in P. sativum, Rai et al. (2016) reported that major QTLs Qruf and Qruf2 were common to all the resistance traits evaluated including the disease severity in the seedling and in the adult plants in different environments.
The minor QTL UpDSIV.2 was environmentally specific and was only detected for DS% in seedlings under controlled conditions. It accounted more than 14% of the phenotypic variation studied. Minor QTLs are observed quite frequently in disease resistance studies, but they are prone to inconsistent expression (Pilet-Nayel et al., 2002;Cobos et al., 2005). The three QTLs (UpDSII, UpDSIV, and UpDSIV.2) showed all negative additive effect. This indicated that the resistance alleles came from the resistant parent, which is supported by the presence of transgressive segregants with a lower disease severity.
Our results have indicated that, taken together, the identified QTLs have explained a very high percentage of the phenotypic variance throughout the population (40.5-51.3% of pea rust severity variation in the field and 62% in controlled conditions), suggesting that an efficient selection could be possible with a few markers tightly linked to the resistant QTLs. The SSR marker AD280 and the DArT-derived marker 3567800 were located at less than 1.5 cM, delimiting a region fairly close to UpDSII. Similarly, the DArT-derived markers 3563695 and 3569323 were located at less than 5 cM of the QTL UpDSIV, as well as the markers 3536422 and 3538798 to the QTL UpDSIV.2. For the QTLs UpDSII and UpDSIV the peak LOD could move 1 cM depending on the year or the conditions of rust evaluation, with some minor changes in the linked markers. This is not unusual, since molecular mapping of disease resistance QTLs in legumes has often revealed the same QTL located in comparable genomic regions but having different markers closely linked depending on how the disease intensity was scored or if the evaluation were carried out in different years (Hamon et al., 2011;Carrillo et al., 2014;Rai et al., 2016). Although the mapping population of our study was evaluated in three different years and different conditions, the genomic positions of the QTLs were comparable in all the cases not without observing that the markers closely linked to the QTLs varied in some of the cases.
Disease resistance genes are commonly organized in complex clusters or loci (Loridon et al., 2005). These regions of the linkage maps are often rich in genes conferring resistance to different pathogens and/or to different races of the same pathogen. After connecting our P. fulvum map with the previously published P. sativum linkage maps (Loridon et al., 2005;Bordat et al., 2011;Carrillo et al., 2014), we have found that the QTLs MpV.2 and MpII.1 conferring resistance to Didymella pinodes (Carrillo et al., 2014) as well as the QTLs Ae-Ps5.2 and Ae-Ps2.2 conferring resistance to Aphanomyces euteiches (Hamon et al., 2011) were co-localized with UpDSII and UpDSIV, respectively. This suggested that the genomic regions of the LGs II and IV of P. fulvum where UpDSII and UpDSIV are likely to control disease resistance and to harbor clusters of disease resistance genes in pea. The identification of such genomic regions involved could be useful in resistance breeding programs through marker assisted selection (MAS).
In this sense, the additional advantage offered by DArTseq TM markers is that they can be readily converted into PCR-based markers (Fiust et al., 2015) in cases where there are still no lowcost markers that closely flank a potential QTL. Also, there is a great possibility of identifying candidate genes that could be used in the foreground selection of favorable alleles since DArTseq produces a large number of markers within gene regions. The identification of potential candidate genes is the first step to identify the genes that control U. pisi resistance in pea. However, additional studies such as functional analysis are needed to finally validate the role of the gene in resistance. The suitability of these genes as candidates for resistance to U. pisi would facilitate an efficient MAS.

CONCLUSIONS
The results obtained from the present study indicated that DArTseq TM provides high-quality markers that can be used to construct dense genetic linkage maps for plants even when no sequence information is available. A total of 37 polymorphic species-specific SSR, STS and SNP markers and 12,021 DArTseq TM based markers were used to develop a reasonably well saturated genetic linkage map of the interspecific P. fulvum RIL population derived from the cross between IFPI3260 × IFPI3251, allowing a precise and fine QTL mapping of important phenotypic traits related to U. pisi resistance in wild pea. The rust resistance loci identified in the LGs II and IV have been a novel report. The expression of the resistance to U. pisi in P. fulvum originated from the resistant parent IFPI3260. DArTseq TM will be very useful both in pea breeding programs and in parallel ongoing projects in legumes.

AUTHOR CONTRIBUTIONS
EB and DR designed the experiments. EB performed the RIL population development and the rust evaluations. EB and MC carried out the P. fulvum DNA extractions as well as the SSRs and STSs evaluation in the RIL population. EC performed the SNPs analysis in the RIL population and critically reviewed the manuscript. EB and MC carried out the QTL analysis and wrote most of the manuscript. DR contributed to the interpretation of results and writing of the manuscript. AK and JC performed DArT analysis, marker selection as well as established the synteny between P. fulvum and M. truncatula. AK and JC also contributed to critical reading.

ACKNOWLEDGMENTS
Authors are greatly indebted to AGL2017-82907-R and AGL2014-52871-R projects co-financed by FEDER for financial support.