Integrated Molecular and Morphological Studies of the Daucus guttatus Complex (Apiaceae)

Abstract 
 In a previous study using 94 nuclear orthologs, we reported the species status of the Daucus guttatus complex to be unresolved, partitioned into three clades. In the present study, a subset of ten of these 94 orthologs was used to infer the phylogeny of the D. guttatus complex and related species. A near parallel set of accessions, planted in a common garden, was used for morphological analyses. The molecular trees are highly resolved for most of the clades, grouping accessions of the D. guttatus complex into four clades. Bayesian concordance analysis and a coalescent approach gave slightly different topologies. Morphological data likewise support four taxa in the complex. Moreover, herbarium research from a companion study informs nomenclature for taxa of the complex. We identify these four clades as D. bicolor, D. conchitae, D. guttatus, and D. setulosus; internested in or among these segregates are the phenetically distinctive species D. glochidiatus, D. involucratus, D. littoralis, and D. pusillus. Our research redefines species variation in the D. guttatus complex, clarifies species names, interspecific relationships, confirms a useful subset of nuclear orthologs for studies of dominant topologies of Daucus, and discovers morphological characters allowing proper identification of the four species of the D. guttatus complex and related species.

Carrot (Daucus carota L. subsp. sativus Hoffm.) is considered among the ten most important vegetables worldwide and is the highest value crop in the Apiaceae (Rubatzky et al. 1999;Simon 2000;Vilela 2004). The genus Daucus is most common in the Mediterranean region and is characterized by dorsally compressed mericarps, hairs on the primary ridges of the mericarp, and spines on the secondary ridges of the mericarps (Okeke 1978).
The latest comprehensive taxonomic monograph of Daucus by Sáenz (1981) recognized 21 species divided into five sections: Daucus sections Anisactis DC., Chrysodaucus Thell., Daucus L., Meoides Lange, and Platyspermum DC. Rubatzky et al. (1999) later estimated 25 species of Daucus. More than 60 species have been proposed as variants within the "D. carota complex" alone (Small 1978), for which there are no or only poorly developed barriers to interbreeding among either the wild forms or their domesticates (Small 1978). Morphological data fail to clearly distinguish many of the taxa within this complex ).
Delineating species boundaries correctly is crucial to understanding the diversity of life (Dayrat 2005). Historically, delineation of species and their taxonomy has been based mostly on morphological characters (Patterson et al. 1993;Pisani et al. 2007). The utility of DNA sequences for taxonomic studies is also well established (Tautz et al. 2003). Different authors have successfully integrated morphological and molecular data to delineate and identify species (Emadzade et al. 2010;Martínez-Flores et al. 2012, Cohen 2014. Morphological and molecular studies have been conducted in Daucus to distinguish taxa with varying results (Vivek and Simon 1999;Bradeen et al. 2002;Spooner et al. 2013Spooner et al. , 2014Arbizu et al. 2014aArbizu et al. , 2014bMezghani et al. 2014;Lee and Park 2014).
Single and low-copy nuclear genes provide significant phylogenetically informative characters complementing other types of data in plant systematics (Zimmer and Wen 2013). Neither nuclear orthologs (Arbizu et al. 2014b) nor morphological data (Spooner et al. 2013) have distinguished the subspecies of D. carota L., but nuclear orthologs and morphology (Arbizu et al. 2014a) generally provide effective data sets to distinguish the other species in Daucus. Relevant to the present study, molecular and morphological studies of multiple accessions of the D. guttatus complex collected in Greece, Iran, Israel, Lebanon, Syria, and Turkey demonstrated high diversity in these taxa, which were partitioned into three wellsupported clades (Arbizu et al. 2014b), and phenetic groups (Arbizu et al. 2014a), suggesting that these putative species needed reevaluation.
Currently, there are no accepted criteria to clearly distinguish the segregates of the D. guttatus complex. They are morphologically similar and difficult to distinguish, causing frequent misidentification of newly acquired germplasm. We here aim to resolve the species boundaries and relationships of these putative species, and other related species as determined by Arbizu et al. (2014b), by integrating molecular and morphological data.

Materials and Methods
Accessions Examined -We examined 83 accessions of Daucus, and one accession of a non-Daucus outgroup Orlaya daucoides (84 accessions in total), collected in 12 countries (Suppl. Data 1). Arbizu et al. (2014b) reported that clade B contains samples of the D. guttatus complex partitioned into three well supported clades (therein called D. guttatus 1, D. guttatus 2 and D. guttatus 3), and other related species that are morphologically distinct from the members of the complex (D. glochidiatus, D. involucratus, D. littoralis, and D. pusillus). Sixty-five of the 84 accessions sampled for the present study correspond to the D. guttatus complex and nine accessions are closely related species (sensu Arbizu et al. 2014b;and including additional sampling). Notably, the epitype of D. guttatus (PI 652342, designated by Martínez-Flores et al. 2016) is included in this sampling. Our morphological analysis included 40 accessions of six Daucus species including most of those examined in the molecular analysis (Suppl. Data 1). All accessions were obtained from the United States National Plant Germplasm System, maintained at the North Central Regional Plant Introduction Station (NCRPIS) in Ames, Iowa. Full details of the accessions examined in this study are available at the Germplasm Resources Information Network -GRIN (http:// www.ars-grin.gov/npgs/acc/acc_queries.html).
DNA Extraction and Sequence Generation -DNA was extracted using the CTAB method (Doyle and Doyle 1990) from young leaves of accessions planted in a greenhouse at the University of Wisconsin-Madison. We used the ten nuclear orthologs that were identified by Arbizu et al. (2014b) as a useful subset to recover dominant topologies in Daucus (the subset determined by the result obtained by the concatenated data set of all 94 orthologs). DNA sequences of 53 accessions of the D. guttatus complex for these ten nuclear orthologs were obtained using Asymmetric PCR Single-Strand Conformation Polymorphism, an efficient alternative technique for isolating allelic variants of Daucus (Rodríguez et al. 2011). Sanger sequencing (Sanger et al. 1977) was conducted at the Biotechnology Center of the University of Wisconsin-Madison. DNA sequences for the other accessions (Suppl. Data 1) were obtained from Arbizu et al. (2014b).
Phylogenetic Analyses -DNA sequences were edited with Staden package version 1.7.0 (Staden 1996) and aligned using MUSCLE version 3.5 (Edgar 2004). Subsequent minor manual alignment corrections and shortening of homopolymers to a maximum of six base pairs (to match the data of Arbizu et al. 2014b) were performed using Mesquite version 3.03 (Maddison and Maddison 2015). We rooted our trees on Orlaya daucoides based on Arbizu et al. (2014b). Maximum parsimony (MP) analysis was conducted in PAUP* version 4.0a145 (Swofford 2002). Missing data and gaps were all scored as missing data. All characters were treated as unordered and weighted equally (Fitch 1971). The most parsimonious trees were found using a heuristic search (Farris 1970) by generating 100,000 random taxon addition sequence replicates using tree-bisection reconnection (TBR) and holding one tree for each replicate. Then, we ran a final heuristic search of the most equally parsimonious trees from this analysis using TBR and MULPARS. Bootstrap values (Felsenstein 1985) for the clades were estimated using 1,000 replicates with a heuristic search, TBR and MULPARS, setting MAXTREES to 1000.
Selection of the best-fit evolutionary models for the individual gene data sets was computed using the Akaike information criterion (AIC), using jModelTest version 2.1.4 (Darriba et al. 2012). With these models, we obtained a maximum likelihood (ML) tree with GARLI version 2.1 (Zwickl 2006). Using the same program, 1,000 non-parametric bootstrap inferences were obtained. Both analyses were performed by the GARLI web service hosted at www.molecularevolution.org (Bazinet et al. 2014).
In addition, we performed a Bayesian concordance analysis (BCA) (Ané et al. 2007) to obtain a primary concordance tree and clade concordance factors with 95% credibility intervals using the program BUCKy version 1.4.3 (Larget et al. 2010). Based on Arbizu et al. (2014b) for the BCA pruned analysis we chose 23 representative accessions (Suppl. Data 1) from major clades as determined from the MP and ML analyses to explore gene to gene conflict in our data set. All ten genes with their corresponding model of nucleotide substitution (Suppl. Data 2) were analyzed separately in MrBayes version 3.2.3 ) using the BEAGLE library  with four chains and two searches run simultaneously for ten million generations sampling every 1,000 generations. This analysis was conducted via the CIPRES (Miller et al. 2010) portal at San Diego Supercomputer Center (http://www.phylo.org). We summarized the MrBayes results for the ten genes using the program mbsum included in BUCKy, removing 1,001 trees from each chain as burn-in. We then performed the BCA with four independent runs with four linked chains for all four different levels of discordance: α = 0.1, 1, 10, and infinite (a larger value of α corresponds to greater gene tree incongruence); in each run with 1,100,000 generations, with 100,000 generations discarded as the burn-in period. Default settings were used for all other parameters.
We also performed a Bayesian analysis using BEAST package version 1.8.2 (Drummond et al. 2012) to obtain a species tree estimation using a coalescent approach. First, an XML format file was generated using BEAUti version 1.8.2 with ten genes with their corresponding model of evolution (Suppl. Data 2) and Yule process as a species tree prior (Gernhard 2008). Based on our MP and ML tree topology, accessions from Daucus clades 1, 2, 3 and 4 (Results) were labeled as different species. Even though our MP and ML results show that D. involucratus is embedded within clade 3, we here considered it as a different species based on the molecular and morphological results of Spalik and Downie (2007), Arbizu et al. (2014a, b), and Lee and Park (2014). Markov chain Monte Carlo (MCMC) chains for both analyses were run for one billion generations sampling every 50,000 generations. We imported the log files into Tracer version 1.6 (Rambaut and Drummond 2014) from the BEAST package to analyze the convergence to the stationary distribution and the effective sample size (ESS) of each parameter. The sample files of trees produced by BEAST were summarized discarding 25% of the trees onto a single "target" tree using the program TreeAnnotator version 1.8.2. The resulting trees were viewed in FigTree version 1.4.2 (http:// tree.bio.ed.ac.uk/software/figtree/).
Morphological Analyses -We considered species that belong to clade B of Arbizu et al. (2014b) for our morphological study. In addition, accessions of Daucus bicolor, not characterized by Arbizu et al. (2014a, b), were considered as part of clade B based on MP and ML analyses of the present study (see Results). Each accession was direct seeded by hand in 1m × 6m observation plots at the West Madison Agricultural Research Station, University of Wisconsin, Madison. We watered the rows by hand and weeded the field with a small plot tiller, hoe, and hand weeding. Low germination rates of some accessions prevented us from having plants in every observation plot.
Daucus accessions were scored and analyzed for 41 characters in total: 31 continuous and 10 nominal characters (Table 1). All characters were always recorded by the same individual and from at least one plant per accession. These characters were selected to represent all those previously used in the latest comprehensive monograph of Daucus (Sáenz 1981) and in subsequent Daucus morphological studies (Small 1978;Arbizu et al. 2014a;Mezghani et al. 2014;Spooner et al. 2014). Plant and leaf characters were evaluated in the field with a ruler or digital calipers, and floral and mericarp characters were measured in the laboratory using a stereoscopic microscope. Forty observation plots lacked plants due to no germination or unintentional elimination during weeding. We obtained one data set including 40 accessions containing all measurements of 82 plants for the 41 characters (plant, leaf, flower and mericarp). Accession PI 652333 (Suppl. Data 1) had neither plant nor fruit data.
We summarized our morphological data using R version 3.2.0 software (R Core Team 2015). For multivariate analyses, means were calculated for the continuous characters using the ddply function in the plyr package (Wickman 2011). A descriptive statistical analysis was conducted to verify the mean, median, standard deviation and range of values. Box plots were used as a graphic tool to visualize comparisons across accessions and to check for outliers that may represent erroneous entries. Accessions were classified into groups based on our molecular phylogenetic analyses. Then we performed stepwise discriminant analyses (linear, common variance) using all 31 continuous variables in our data set comprising 40 accessions.
Discriminant analyses were performed in JMP ® version 11.2 (SAS Institute Inc. 1989-2007) using a backward selection; a model with significant variables in identifying accession structure was obtained by removing characters one at a time until the model F-test p value ≤ 0.05. We then performed canonical variate analyses (CVA). Hierarchical cluster analysis (HCA) was conducted using all 10 nominal characters and the significant continuous characters previously obtained. For our data set (31 nominal characters and 40 accessions in total), we considered the following: (1) all taxa with 31 nominal characters (plant, leaf, flower and mericarp); and (2) all taxa with a subset of charactersthe mericarp characters only. Analysis using the mericarp characters only was performed to see if mericarps alone could be used to distinguish the species analyzed in this study. In addition, a linear discriminant analysis (LDA) was performed with the data mentioned above using the lda function in the MASS package (Venables and Ripley 2002) in R. Results were very similar to our CVA and we do not discuss this further.
In concert with our companion study (Martínez-Flores et al. 2016; this issue of Systematic Botany), we studied the original descriptions and the floristic literature to assign species names to the four molecular clades (putative species) of the D. guttatus complex. Proper names were assigned by (i) analyzing digital pictures from GRIN website and our experimental plots, and (ii) checking names with herbarium specimens. Potential synonyms were inferred through: (1) the latest comprehensive taxonomic treatment of Daucus (Sáenz 1981); (2) regional floras outlined in the Introduction; and (3)

Results
Completeness of Sequencing -Our aligned DNA data set consisted of 10,982 base pairs (bp). The alignment, Supplementary Data 1 and 2, and Supplementary Figs. 1-4 are deposited in Dryad (http://datadryad.org/review?doi=doi:10.5061/ dryad.mp10n). Sequences of the 53 accessions newly sequenced here (see Suppl. Data 2 for GenBank accession numbers) resulted in one or two alleles differing only in minor differences in every gene. In previous work, the topologies of the phylogeny of Daucus using a single allele vs. two alleles merged into one showed minor topological differences (Arbizu et al. 2014b). As a result, we analyzed the one allele containing the least number of missing data (base pairs). Of these ten markers/53 accessions matrices (530 sequences) there were 48 missing sequences, resulting in 9.1% missing data and we did not observe any ambiguous alignments. The entire data matrix of 840 sequences (ten markers/84 accessions matrices) had 57 missing sequences (6.8% missing data).
Maximum Parsimony Analysis -Concordant with prior molecular studies in Daucus (Spalik and Downie 2007;  Figure 1). Clade A' contains all species possessing 2n = 18 chromosomes, while the remaining taxa in clades A and B possess 2n = 16, 20, 22 and 44 chromosomes, and the cladistic relationships match those of Arbizu et al. (2014b). Within clade B, our MP tree recovered four major clades, three of which were found by Arbizu et al. (2014b), and one additional clade including samples not previously studied. All four of these D. guttatus species complex clades had high bootstrap values of 98-100%. Like the results of Arbizu et al. (2014b), all accessions of D. involucratus form a clade, but in the present analysis they are internal to other members of clade 3 rather than sister to it. Also unlike Arbizu et al. (2014b), the single accession of D. glochidiatus is sister to D. littoralis and other members of D. guttatus clade 1 in the present study, whereas D. glochidiatus was sister to D. pusillus and other members of D. guttatus clade 2 in Arbizu et al. (2014b). All other species relationships between this study and Arbizu et al. (2014b) are the same.
The use of additional accessions of the D. guttatus complex in this study revealed a geographical component with the two accessions of D. guttatus clade 1 collected in Israel (PI 279763 and PI 295858; 100% bootstrap). Also, clade 2 contains all accessions collected in Greece, except for two accessions collected in Turkey (Ames 25814, PI 652360). Clade 3 contains accessions only collected in Turkey, except PI 652332 (D. involucratus).
Maximum Likelihood Analysis -The ML tree using separate models for every single gene with 1,000 bootstrap replicates ( Fig. 1) has the same overall topology as the MP tree (Suppl. Figure 1), including the geographic subset in clade 1. It also resolves the same clades, A, A' and B. Furthermore, there are good bootstrap support values in most branches of this tree. As in the MP analysis, D. involucratus was embedded in clade 3.
Bayesian Concordance Analysis -Results from our pruned analysis showed good mixing and good results as the standard deviation of concordance factors was less than 0.005. The primary concordance tree (Fig. 2), estimated for ten genes and 23 accessions with Bayesian concordance analyses, showed two differences in the topology compared to the MP and ML trees: (1) D. glochidiatus is sister to a clade formed by accessions of clades 1, 3, 4, and D. littoralis; (2) D. pusillus is sister to all species mentioned in (1). Daucus pusillus and D. glochidiatus are supported by one gene as sister clades of the species mentioned above. Daucus involucratus was also embedded in clade 3. There were no differences among the concordance factors using the four different prior probabilities on gene tree incongruence (α value). Clades A and B are supported by 7.2 and 7 genes, respectively. Less than four genes support clades 1, 2, and 4. There are 5.4 genes supporting clade 3.
Species Tree Estimation -Under the multispecies coalescent approach, the Bayesian analysis using 84 accessions containing 18 species and ten genes with their corresponding model of evolution gave very high ESS values (higher than 1000) for most of the parameters. Concordant with MP, ML, and BCA tree topologies, there are two main clades of Daucus, A and B (Fig. 3). Values of support (i.e. posterior probabilities) are slightly higher for many branches (Fig. 3) compared to the species tree of Arbizu et al. (2014b). The species tree resulted in a similar topology to the MP and ML. However, there are three notable exceptions: (1) D. glochidiatus is no longer sister to clade 1 and D. littoralis, but rather, with low support, to clade 3 (D. conchitae with D. involucratus); (2) Clade 2 and D. pusillus form a clade with low posterior probability with clade 1, D. littoralis, and clade 4; and (3) D. aureus and D. muricatus form a single clade supported by a high value of support.
Morphological Analyses -Our morphological analyses considered the accessions grouped in the four discussed clades within clade B as distinct taxa as determined by the MP and ML analyses (except that D. involucratus, as discussed, was considered separately). Graphical analyses of the 41 character state distributions examined in this study are shown in Suppl. Figure 2. Patterns of variation were very similar to the studies reported by Spooner et al. (2014) and Arbizu et al. (2014a). That is, while a few characters showed little to no variation, most characters showed tremendous variation within all taxa and overlap of ranges across species. However, a combination of morphological characters was useful in identifying species within clade B.
Twenty-two of the 31 continuous variables were significant in the F-test, p ≤ 0.05 in the stepwise discriminant analyses in at least one of the two analyses. The following nine were not significant: leaf width at maximum, pigmented central umbellules diameter, bract length, bract width, length of shortest peripheral ray, peripheral petal length, stamen length, stylopodium length, and vittae width.
Canonical variate analysis was reported to be a useful discriminant method to distinguish taxa of Daucus Arbizu et al. 2014a). Furthermore, Arbizu et al. (2014a) showed that HCA produced a tree based on average similarity of all characters of the three subsets of D. guttatus. Two different analyses were conducted using discriminant approaches and we present the F-test p values of the characters retained in all stepwise discriminant analyses (Table 1). Canonical variate analysis separated all taxa easily into six groups: the four D. guttatus complex clades, D. involucratus and D. littoralis (Fig. 4A). Stepwise discriminant analysis identified 21 of the 31 continuous characters as significant discriminators for this data set within all taxa at the F-test p value ≤ 0.05 (Table 1, column 3). Character state distributions of the six most significant discriminators are displayed in Fig. 5 and as expected, character state overlap even within these best discriminators was common. Eight of the 21 significant characters are mericarp characters. However, an analysis of only mericarp characters failed to clearly separate them into their corresponding groups (Fig. 4B), except accessions from clade 3, D. involucratus and D. littoralis.
Stepwise discriminant analysis identified five of the 10 continuous characters as significant discriminators within all taxa at the F-test p value ≤ 0.05 (Table 1, column 4). Only one accession of clade 4 (Ames 25833) is intermingled within clade 2. Separation of accessions of clades 1, and 2 is not obvious since there is one accession of these two clades placed in other groups (Fig. 4B). The HCA shows a different grouping pattern. It did not cluster all taxa with clear resolution. Daucus littoralis is considered the most distinctive species (Suppl. Figure 3).
Species Names in the Daucus guttatus Complex -Based on a comparison of our molecular and morphological results with the original species descriptions, and the 482 SYSTEMATIC BOTANY [Volume 41    Primary concordance tree of the Daucus guttatus complex including species that belong to Clades A and B obtained with Bayesian concordance analysis using 10 nuclear orthologs and 23 accessions. Numbers above the branches are the concordance factors, which do not show significant differences for different α values (0.1, 1, 10, and infinite). Clades 1, 2, and 3 were identified in Arbizu et al. (2014b). Clade 4 is a new finding in the present study.     Arbizu et al. (2014b) demonstrated that a subset of ten of 94 nuclear orthologs provided effective discriminators for species boundaries and relationships of Daucus but highlighted unresolved species boundaries in the D. guttatus complex that were partitioned into three subclades within clade B. Our present study explored this problem through the use of additional accessions of the D. guttatus complex and discovered them to now be partitioned into the original three clades and one additional clade (now four). Morphological analyses of these accessions provided support for these as four species of the D. guttatus complex (D. guttatus, D. bicolor, D. setulosus and D. conchitae), but with characterstate overlap among the species. Concatenated analyses and BCA placed all accessions of D. involucratus together, but internested entirely within D. conchitae (clade 3 of the D. guttatus complex). This can be explained by the use of a reduced number of nuclear orthologs (10 vs. 94) compared to the study of Arbizu et al. (2014b).
The primary concordance tree (PCT, Fig. 2) produced a topology slightly different from MP and ML analyses. BCA does not assume that genes have the same topology (Ané et al. 2007). Concordance factors (i.e. genomic support) are less than 50% for the clades 1, 2, and 4, showing discordant gene histories in the genomes of Daucus, as reported by Arbizu et al. (2014b). On the other hand, the PCT analysis shows very high concordance factors for the clade containing D. carota subsp. capillifolius and D. sahariensis. In addition, the clade grouping two accessions of clade 1 (PI 652387 and PI 652343), and clade 2 (Ames 25636 and PI 652325) also have high concordance factors indicating less discordant histories in these genomes. The most likely factors that may be causing discordance are recombination, hybridization and introgression (Rieseberg et al. 2000), and incomplete lineage sorting (Pamilo and Nei 1998). The species tree topology is similar to the PCT.
The multispecies coalescent approach (Fig. 3) assumes that genealogical discordance is entirely due to incomplete lineage sorting, the most common cause of discordance among gene trees (Maddison and Knowles 2006;Degnan and Rosenberg 2009;Baum and Smith 2013). We place greater weight on the BCA results since it makes no assumption about the reason of genealogical discordance (Larget et al. 2010). The species tree obtained by Arbizu et al. (2014b) included 104 accessions of 27 species of Daucus and outgroups and showed a topology that was not congruent with the concatenated analyses of that study. In the present study, we report a species tree different from the species tree of Arbizu et al. (2014b). This can be explained by the following factors: (1) ESS values for all parameters are very high in our present study (higher than 500) indicating good mixing; and (2) we used ten selected genes that resemble the dominant topology of Daucus.
Morphological Analyses -Discriminant analysis was used by Spooner et al. (2014) and Arbizu et al. (2014a) to explore the validity of species and subspecies within Daucus, showing good discrimination except for the subspecies of 488 SYSTEMATIC BOTANY [Volume 41 D. carota. The taxonomic value of fruits in the Apiales have been reported by Lee et al. (2001), Liu (2004), Baranski et al. (2006), Liu et al. (2006) Akalin Uruşak and Kizilarslan (2013), Mezghani et al. (2014), and Tavares et al. (2014). Sáenz (1981), Pujadas Salvà (2003), Baranski et al. (2006), and Mezghani et al. (2014) indicated that fruits provided critical taxonomic characters distinguishing taxa within Daucus since fruits tend to show less phenotypic plasticity in comparison with other plant organs. Other studies questioned the utility of fruit characters alone as a practical tool to classify taxa in the Apiales Katz-Downie et al. 1999;Spalik et al. 2001;Vallejo-Roman et al. 2002;Lowry et al. 2004). We investigated this issue in the present study by examining only the 10 continuously variable fruit characters. Our discriminant analyses identified a subset of five mericarp characters that are most useful to distinguish D. conchitae, D. involucratus and D. littoralis (Table 1, Fig. 4B), but not the other species. Thus, a combination of mericarp and other vegetative and floral characters are needed to distinguish members of the D. guttatus complex (but with some character overlap). Using 11 accessions from the D. guttatus complex, Arbizu et al. (2014a) could separate three species of this complex using HCA. On the other hand, in the present study our HCA could not separate the four species of the Daucus guttatus complex (Suppl. Figure 3). This difference can be explained by the greater number of accessions (37 in total for the D. guttatus complex) used in the present study.
Association of Species Names to the Daucus guttatus Complex (clades 1, 2, 3, and 4) -To date, the species boundaries and application of existing names to members of the D. guttatus complex have remained unclear. Our molecular analyses identified four well-supported clades of members of this complex and showed that they had morphological support, but only with multiple characters that exhibited some overlapping sets of character states (Suppl. Figure 2). Members of clade B possess chromosome numbers of 2n = 20, 22 and 44, but many of them are from few counts. The accuracy of 2n = 20 and 22 counts and their possible association with clades in the D. guttatus complex need further investigation with additional accessions.
The original descriptions of D. bicolor and D. guttatus by Smith (1819) discriminated these two species by differences in plant height, degree of branching, pubescence, bract size, and color pattern of central umbellets. Daucus broteri was originally described by Tenore (1830) as a plant with hairy leaves and tripinnate, compact umbels, multiple umbel rays, almost all bracts 3-fid, wide fruit spines confluent at the base, and 6-8 spines per ridge. Greuter (1979) described D. conchitae based on plants without fruits. He described it as a plant branched at the base, with five to seven simple bracts with a length of 15 mm, exceeding the umbel perimeter. He mentioned a distinguishing character of D. conchitae as simple or entire bracts. Daucus setulosus (de Candolle 1830) was described as a plant with stems branched and setulose, leaves pinnatisect with various segments, and bracts with multiple segments. Different authors have considered D. setulosus as a synonym of D. guttatus and/or D. bicolor. We here quantify tremendous variation in all of the above characters. Our analysis showed that these species are similar, with sometimes overlapping character states, but with the best discriminating characters as shown in our key to be involucral bract length and shape, primary umbel diameter, coloring of the central umbellets, mericarp size, vittae shape, and degree of confluency of the spines. Clearly, discrimination of species in clade B is difficult and relies on the use of a suite of characters that overlap in range (a polythetic approach). Ultimately, we associate species names with the clades by examining characters of the majority of accessions in these clades and comparing to type specimens (Martínez-Flores et al. 2016).
Flower color variation is a model for the experimental study of evolution. Clegg and Durbin (2000) found that almost all mutations that determine phenotypic differences in Ipomoea purpurea (L.) Roth flower color variations are the result of transposon insertions. This may be a plausible explanation for the case of dark-colored pattern present in D. guttatus clades 1, 2, 4. Westmoreland and Muntan (1996), and Goulson et al. (2009) favored the hypothesis that the dark flowers in umbels of Daucus significantly increase both insect attraction and fruit production. Conversely, Polte and Reinhold (2013) demonstrated that the dark central floret in wild Daucus has no role in pollination, but may play a role in reducing parasite infestation by the gall midge Kiefferia pericarpiicola (Bremi, 1847). This issue needs further research to clarify the role of purple-colored flowers in Daucus species.
Most of the accessions examined in clade 1 were collected in Syria very close to the collection conducted by Al-Safadi (2008) (Slanfah, Salma and Rabiaa, and near Tartus). Interestingly, accessions that matched the original description of D. guttatus were placed in clade 1, including the epitype of D. guttatus (PI 652342). Accessions of clade 4 match the original description of D. bicolor. They look similar to D. bicolor due to the dark-colored pattern in the umbels, but bracts do not always exceed the umbel perimeter. However, bracts are 3-fid. Peripheral petal length is longer than in D. guttatus, D. setulosus and D. conchitae. Leaf, petiole, and peduncle pubescence is smoother. Mericarp length and width is shorter and has a triangular vittae shape. As a consequence, we identified clades 1 and 4 as D. guttatus and D. bicolor, respectively.
There are accessions in Italy that have been reported as D. broteri and were originally misidentified as D. muricatus because of the wide mericarp spine confluent at the base (Martínez-Flores et al. 2016). Based on morphological evidence only, it is possible that those accessions reported in Italy may be D. guttatus or yet another species in this complex. Sáenz (1981) classified D. muricatus and D. bicolor, and D. guttatus within sections Platyspermum and Daucus, respectively. She considered D. broteri and D. setulosus as synonyms of D. bicolor and D. guttatus, respectively. We suggest that samples examined by Sáenz (1981) for D. bicolor are one form of D. guttatus with spines strongly widened at the base, confluent into a crest (the same as Ames 25733, PI 652340, among others). In addition, Sáenz (1981) described D. guttatus as a species with bracts divided and exceeding the umbel perimeter, spines slightly confluent at the base, and large triangular shaped vittae. However, we propose samples examined are D. setulosus as they match the bract description reported by de Candolle (1830) and our mericarp description for clade 2. Therefore, we conclude accessions of clade 2 correspond to the species D. setulosus.
Daucus involucratus is embedded within clade 3 and was previously reported as a close relative to D. guttatus 3 2016] ARBIZU ET AL.: STUDIES OF THE DAUCUS GUTTATUS COMPLEX (Arbizu et al. 2014b; here named as clade 3; chiefly D. conchitae). Greuter (1979) indicated that D. conchitae resembles D. involucratus in general habit. On the contrary, Sáenz (1981) indicated D. conchitae appeared to be a taxon close to D. littoralis because of their anatomical characters. Aytaç and Duman (2013) collected individuals in Antalya, Turkey, with simple bracts and identified them as D. conchitae. Constantinidis (2013) evaluated fruit morphology of all plants collected in the Kastellorizo Island group and concluded D. conchitae is simply a form of D. guttatus. He reported his collection has individuals with simple and pinnatifid bracts, usually in the same umbel, and fruit spines weakly confluent at the base. His descriptions are very similar to the accessions examined here placed in clade 3. Surprisingly, molecular analyses placed D. conchitae as a close relative of D. involucratus (Spalik and Downie 2007;Lee and Park 2014). We here infer that most accessions of clade 3 belong to D. conchitae (with specimens we identify as D. involucratus also nested within this clade, as discussed). Vittae shape has been considered as a discriminant character within taxa of Daucus (Thellung 1926;Sáenz 1981;Heywood 1982). We documented the mericarp of D. conchitae to have elliptical vittae. This characteristic was not present in any other species examined here, except the accessions of D. guttatus with large spines confluent at the base forming a wide crest, and D. littoralis. The use of wild crop relatives in plant breeding is likely to intensify as breeding techniques improve and adaptation to climate change becomes more pressing (Khoury et al. 2013). There is a constant demand for development of new hybrid carrot cultivars (Rubatzky et al. 1999;Frese and Nothnagel 2008) and wild Daucus species have utility for breeding against pest and disease tolerance or resistance, yield increase, male sterility, and for nutraceutical and culinary traits (Camadro et al. 2008). To better use germplasm maintained at different germplasm banks, it is imperative to determine the species boundaries and phylogenetic relationships of Daucus. Wild Daucus species possessing 2n = 18 chromosomes are not the only source of genetic resources for carrot. Daucus pusillus Michx. (2n = 22), a wild carrot species from the American continent, was tested to establish the feasibility of hybridization with the cultivated carrot under controlled pollination, resulting in the possibility of obtaining interspecific hybrids (Camadro et al. 2008). Daucus guttatus, with 2n = 20 chromosomes, may also be used as a potential breeding source. It has been reported to have prominent antibacterial activity against Corybacterium pyogenes (Radulovíc et al. 2011), and (as D. broteri) to have resistance to carrot fly, Psila rosae (Fabricus, 1794) (Hardman et al. 1990). Daucus littoralis (2n = 20), a species closely related to D. guttatus, shows anti-mycobacterial activity against one of the most lethal infections worldwide, tuberculosis (Başer et al. 2008). Indeed, this species can be considered a good candidate for the potential production of new carrot hybrids with new useful bioactive components.