Twenty eight carrot accessions comprising four wild carrots of different origin, 23 western type carrot cultivars representing four types of root shape and a DH1 plant (Iorizzo et al., 2016) as the reference, were used for ILP validation (Table 1). Total genomic DNA was isolated from fresh young leaves using commercial DNeasy Plant Mini Kit (Qiagen) and used as the template for PCR amplification.

Coordinates of 4028 DcSto insertions belonging to 14 families were compared to coordinates of ca. 32 thousand genes annotated in the carrot reference DH1 genome assembly (Iorizzo et al., 2016; NCBI accession LNRQ01000000). 609 gene-associated DcSto insertion sites localized in introns were identified, of which 209 were manually selected for development of ILP markers. The criteria for initial selection were as followed: insertion sites were (1) free from any other annotated repetitive sequences, (2) present in introns not longer than 3.7 Kb, and (3) evenly distributed over each chromosome. Primer3 (Untergasser et al., 2012) and Primer-BLAST (Ye et al., 2012) were used to design PCR primer pairs anchored in exons flanking introns harboring the selected DcSto insertions. Primer pairs were designed to amplify fragments in a 400–3,700-bp range. The optimal annealing temperature was set to 58°C; and the size and GC content ranged from 18 to 23 bases and 40 to 60%, respectively.

Candidate ILP markers were selected for experimental evaluation. Amplification was carried out in a 10 μL total volume containing 20 ng of genomic DNA, 0.5 μM each of forward and reverse primer, 0.25 mM of each dNTP (Thermo Fisher Scientific), 0.5 U Taq DNA polymerase (Thermo Fisher Scientific) and 1x Taq buffer. The PCR amplifications were performed in an Eppendorf MasterCycler Gradient using the following thermal profile: 94°C (120 s), 30 cycles of 94°C (30 s), 56°C (30 s), 68°C (120 s) and final step of 68°C (600 s). For primers generating ambiguous profiles, the annealing temperature was adjusted to 58, 59, or 60°C. PCR products were separated in 1% agarose gels run in 1x Tris-borate-EDTA buffer (pH 8.0) at a constant current of 5V/cm for about 2 h, stained with Midori Green (Nippon Genetics) and analyzed using GelDoc-It imaging system (UVP). GeneRuler 1 kb and 100 bp⁺ DNA Ladders (Thermo Fisher Scientific) were used to determine product sizes for each locus. The amplicons representing additional local rearrangements within introns were excised, purified using GenJET™ Gel Extraction Kit (Thermo Fisher Scientific), cloned into T/A cloning vector (Promega Corporation) and transformed into Escherichia coli, strain DH10B. Up to five recombinant colonies were selected and cultured overnight at 37°C in culture tubes containing 5 mL of Luria–Bertani medium and ampicillin (100 mg/L). Plasmids were purified using Wizard SV Minipreps KIT (Promega Corporation). Sequencing reactions were set up with universal primers sp6 and T7 using Big Dye terminator chemistry (Applied Biosystems), as recommended by manufacturer. Sequencing was carried out on ABI 3700 capillary sequencer (Applied Biosystems).The sequences were manually edited using BioEdit (Hall, 1999) and aligned to the sequences of predicted genes for which ILP primers were designed.

The ILP marker profiles were scored manually. Each allele was scored as: 1 (empty insertion site), 2 (occupied insertion site) or 0 (lack of amplification).The codominant marker matrix with diploid individuals was created (Supplementary Table 1) and used in GenAlEx 6.5 (Peakall and Smouse, 2006) for creating genetic distance matrix and analysis of molecular variance (AMOVA). Expected and observed heterozygosity (H_e and H_o), and fixation index (F_IS) were computed using POPGENE 1.32 (Yeh et al., 2000). Polymorphism informative content (PIC) of n-allele locus, an indicator of a genetic marker's usefulness introduced by Botstein et al. (1980), was calculated as: PIC=1−∑i=1npi2​−∑i = 1n−1∑j = i+1n2pi2pj2​, where p_i and p_j are the population frequency of the ith and jth allele. Genetic structure was inferred using Bayesian model-based software STRUCTURE 2.2.3 (Pritchard et al., 2008) without information on the accession origin. Ten independent iterations with an admixture and correlated allele frequencies model were performed. The length of the burn-in period and the number of Markov Chain Monte Carlo (MCMC) replications after the burn-in were assigned at 10⁵ for each number of clusters (K) set from 1 to 27 and 1 to 23 for further subclustering. The estimation of K was provided by joining the log probability of data [LnP(D)] from STRUCTURE output and an ad hoc statistics ΔK (Evanno et al., 2005) based on the second rate of change of the log probability of data with respect to the number of clusters. In addition, CLUMPAK software (Kopelman et al., 2015) was used to confirm the selection of the best K. Based on the chosen K, each carrot accession was assigned to a subpopulation for which its membership value (Q) was higher than 0.6. AMOVA was performed using GenAlEx 6.5 to evaluate differentiation among the subpopulations. Principal coordinate analysis (PCoA) was conducted to visualize genetic diversity of the studied accessions.

Insertion sites of 209 DcSto MITEs within introns of annotated genes were chosen to develop Daucus carota Stowaway-like Intron Length Polymorphism (DcS-ILP) markers evenly distributed throughout the genome (Figure 1). The number of DcSto insertion sites evaluated per chromosome varied from 18 (chromosome 9) to 32 (chromosome 2), with an average of 23.22. Their density ranged from 1.37 (chromosome 2) to 2.57 per Mb (chromosome 1), with an average of 1.76.

Upon PCR amplification, 100 of the 209 sites showed the expected DcSto insertion-based polymorphism, however, in case of 10 sites at least one additional amplicon was present in at least one accession (Figure 2). Sequencing of those amplicons revealed that none of the additional variants was related to the activity of the DcSto copy present in the reference genome (data not shown). Of the remaining 109 sites, six did not amplify efficiently; 32 were monomorphic for all tested plants; 13 showed a complex pattern resulting from nonspecific amplification, whereas 58 yielded polymorphic products not associated with DcSto insertions (i.e., sizes of PCR products did not correspond to the expected sizes of empty or occupied variants) (Table 2).

The length of introns harboring the selected DcSto insertions varied from 449 to 3,637 bp. Based on the length of amplified introns, the developed markers were divided into six classes; I to V with intron size ranging from 400 to 3,400 bp, each at 600-bp interval, and class VI comprising introns longer than 3,400 bp (Table 3). Introns belonging to classes I to IV comprised 97.6% of all the developed markers. Class I and II markers were the most numerous, whereas class III markers showed the highest (55.6%) successful amplification rate indicating the most suitable length of introns considered for ILP markers. DcS-ILP markers of class V and VI were characterized by ambiguous amplification patterns, therefore not considered for further analyses.

Finally, 90 DcS-ILP (Supplementary Table 2) markers showing biallelic DcSto insertion polymorphism (Figure 2) were chosen for development of a panel for genotyping the carrot.

The utility of 90 biallelic DcS-ILP markers was verified by estimating the genetic diversity of the collection of 27 D. carota accessions comprising 23 cultivated and 4 wild populations. In total, 180 alleles were identified with an average of 2.0 per locus. 2.78% of the alleles were rare (frequency <0.05) and the mean effective number of alleles was 1.56. The observed heterozygosity for individual loci ranged from 0.04 to 0.56, with an average of 0.24, whereas the expected heterozygosity ranged from 0.04 to 0.51, with an average of 0.34. Shannon's index was from 0.09 to 0.69, with an average of 0.50. Among all the loci analyzed with the Wright's fixation index, 67 were positive. The PIC values ranged from 0.04 to 0.37, with an average of 0.27 (Supplementary Table 1).

STRUCTURE analysis based on 90 loci representing DcSto insertion-derived polymorphisms was performed to evaluate genetic structure of the 27 accessions. The value of ΔK statistics was the highest when two clusters were assumed [ΔK₍₂₎ = 297.64]. The increase in the number of assumed clusters resulted in low ΔK value [ΔK_(>2) = 0.01–52.35]. Twenty three cultivated accessions were assigned to cluster 1 (C1) with membership coefficients (Q) ranging between 0.831 and 0.997, whereas cluster 2 (C2) comprised exclusively wild accessions with the Q value of 0.965–0.998 (Figure 3A). The level of genetic diversity within C1 (0.31) was slightly higher than within C2 (0.29).

To evaluate the genetic structure of the 23 cultivated accessions further subclustering was performed on the accessions assigned to C1. The highest ΔK was observed for K = 21 [ΔK₍₂₁₎ = 22.77], K = 2 [ΔK₍₂₎ = 17.33] and K = 4 [ΔK₍₄₎ = 14.55]. ΔK values for K = 3, K = 5–20 and K = 22–23 were not significant (ΔK = 0.164–4.16). The mean value of log probability of the data was higher for K = 4 than for K = 21, and K = 2 [LnP(D)_{K = 4} = −1891.7, LnP(D)_{K = 2} = −1922.5, LnP(D)_{K = 21} = −2703.2], therefore four subclusters were chosen as the most probable genetic structure of the studied cultivated accessions. With K = 4, three accessions were assigned to subcluster SC1 with Q ranging between 0.928 and 0.962, six to subcluster SC2 with Q between 0.746 and 0.908, five to subcluster SC3 with Q between 0.825 and 0.954 and five to subcluster SC4 with Q between 0.782 and 0.922 (Figure 3B). Four accessions, namely Chantenay Red Cored, Chentenay Rex RS, Danvers 126, and Danvers could not be assigned to any of the subclusters due to high level of admixture (Q < 0.6). The overall Q proportion of each of the four types clearly distinguished (Q > 0.6) the membership of Chantenay root type in SC1 (Q = 0.605), Danvers root type in SC2 (Q = 0.626), Imperator root type in SC3 (Q = 0.785), and Paris Market root type in SC4 (Q = 0.884) (Table 4).

AMOVA attributed 19% (P = 0.001) of the total genetic diversity to variation among the root types. The diversity of the 23 cultivated accessions was revealed by PCoA (Figure 4). Using the first three axes 31.7% of the total variation could be explained, with the 1st, 2nd, and 3rd axes explaining 12.1, 10.4, and 9.2%, respectively.

The above results suggested four separate groups in the collection of 23 cultivated carrots and the grouping generally corresponded with a postulated breeding history of western carrot types presented by Banga (1963), indicating that Chantenay and Danvers types originated from the Late Half Long Horn group, while Paris Market type descended from the Early Short Horn group. Both historical groups differ in terms of their storage root shape and earliness. In turn, the origin Imperator type was traced back to a cross between Chantenay and Nantes (Figure 3C).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.