Taxonomy of the traditional medicinal plant genus Ferula (Apiaceae) is confounded by incongruence between nuclear rDNA and plastid DNA

The aim of this study was to identify major clades in the economically important and taxonomically difficult genus Ferula (Apiaceae tribe Scandiceae subtribe Ferulinae) to provide a classification framework. Phylogenetic relationships among 126 of Ferula spp. and eight species of its sister genus Leutea were evaluated based on nuclear ribosomal (nr) DNA ITS and three plastid regions: the rps16 intron, the rpoC1 intron and the rpoB-trnC intergenic spacer. One hundred and fifty-three accessions were considered including type specimens of seven species. Congruence between nrDNA and plastid DNA data was assessed using a hierarchical likelihood ratio test. Phylogenetic trees were inferred using maximum likelihood and Bayesian methods. Terminals introducing topological conflict were ascertained using two approaches: identifying (1) these with significantly different positions between nrDNA and plastid DNA trees; and (2) a set of rogue taxa in combined trees that, when removed, increased tree resolution and bootstrap support. The results demonstrate significant incongruence between nrDNA and plastid DNA data that persisted after the removal of 41 terminals identified as having significantly different position in nrDNA and plastid DNA trees or after the removal of 13 terminals identified as rogue taxa. Comparison of nrDNA and plastid DNA trees suggest intensive reticulate evolution, particularly in the Irano–Turanian floristic region. Traditional classification systems of the genus are not supported by molecular data, whereas some lineages apparent in molecular trees, particularly Chinese and Mediterranean endemics, are congruent with morphological characters and/or biogeography. We select lectotypes for several infrageneric names and propose a new classification system of the genus with four subgenera and ten sections.


INTRODUCTION
Among the genera of Apiaceae, Ferula L. is one of the largest, includes many species commonly used in traditional medicine and is a promising source of biologically active ingredients. However, the current classification system of the genus does not match the phylogenetic trees, providing little support for research and commercial use. With many species, a wide distribution and poorly understood morphology, the genus epitomizes all the nightmares of a taxonomist. The c. 170 species of this Eurasian *Corresponding author. E-mail: spalik@biol.uw.edu.pl genus are distributed from the Canary Islands in the west through the Mediterranean region, Middle East and Central Asia to western China in the east and northern India in the south (Korovin, 1947;Pimenov & Leonov, 1993). The genus has not been taxonomically revised since the monograph of Korovin (1947). Regional floristic treatments resolved some critical issues but also contributed to the taxonomic confusion because they were inconsistent with one another (Korovin, 1951(Korovin, , 1959Peşmen, 1972;Safina & Pimenov, 1984;Chamberlain & Rechinger, 1987;Safina & Pimenov, 1990). Substantial morphological variation in the genus has resulted in the description of multiple nearly indistinguishable species, which are sometimes represented in herbarium collections by few, poorly preserved specimens. Fruits and basal leaves are essential for correct identification (Peşmen, 1972); the latter, however, are often absent at fruit maturity in monocarpic taxa. Because of the large size of these plants, herbarium specimens usually include only lateral branches and lateral divisions of basal leaves. An experienced taxonomist may identify most species from flowering material (Chamberlain & Rechinger, 1987), but, as underlined by Korovin (1947), the circumscription of species should be based on observations of living plants because it requires examination of complete specimens with roots, stem bases, basal leaves, inflorescence, flowers and ripe fruits. Given these constraints, a taxonomic revision of Ferula seems currently impossible for a single taxonomist. Therefore, it would be useful to subdivide the genus into smaller, monophyletic and workable units that could be subject to partial revisions.
In his monograph, Korovin (1947) established several subgenera and sections based mostly on habit and vegetative features. Safina & Pimenov (1983,1984,1990 contested these infrageneric divisions and proposed an alternative classification system inferred from fruit morphology and anatomy. However, subsequent molecular studies have not supported either of these treatments, indicating that the infrageneric classification system of Ferula has to be constructed anew (Kurzyna-Młynik et al., 2008;Panahi et al., 2015).
Ferula was traditionally classified in tribe Peucedaneae (Pimenov & Leonov, 1993). However, immunochemical studies demonstrated that the genus is not closely related to other genera in this tribe (Shneyer, Borschtschenko & Pimenov, 1995;Shneyer, Kutyavina & Pimenov, 2003). Phylogenetic studies using nuclear ribosomal internal transcribed spacer (nrDNA ITS) sequence variation revealed that Ferula is placed in tribe Scandiceae and forms a clade with Dorema D. Don and Leutea Pimenov (Kurzyna-Młynik et al., 2008). This clade was formally recognized as subtribe Ferulinae. However, the resolution of the ITS tree was low and the infrageneric subclades received moderate internal support. Subsequent analyses using ITS and three non-coding plastid DNA (plastid DNA) sequences showed that Leutea is sister to Ferula, whereas Dorema is nested in it; consequently, Dorema was subsumed under Ferula (Panahi et al., 2015). However, the phylogenetic signal differed significantly between nrDNA and plastid DNA markers resulting in a poor resolution of the tree at the infrageneric level. Moreover, this study included only 68 Ferula spp. of > 170 recognized in the genus.
Here, we further investigate the phylogenetic relationships in Ferulinae using a comprehensive sampling of species and analysing data from nuclear and plastid genomes. Our aim is the identification of major clades of Ferula that may be subject to partial taxonomic revisions and, therefore, provide a phylogenetically meaningful classification system of the genus.

laboraTory procedures
Total genomic DNA was extracted from c. 10-20 mg of dried plant material using the DNeasy Plant Mini Kit (Qiagen, Valencia, CA, USA) following the manufacturer's instructions. Internal transcribed spacer regions (nrDNA ITS) and three non-coding plastid DNA regions (the rps16 intron, the rpoC1 intron and the rpoB-trnC intergenic spacer) were chosen for phylogenetic reconstruction. These markers have already been used in our previous analyses of Ferulinae (Panahi et al., 2015).
PCR products were separated using 1% agarose gel with ethidium bromide or Midori Green DNA Stain (ABO) and purified using the QIAEX II Agarose Gel Extraction Kit (Qiagen). Sequencing from both strands was performed using Big Dye terminators (Applied Biosystems, Foster City, CA, USA). The chromatographs were assembled and edited using SeqMan Pro ver. 12 (Dnastar, Madison, WI, USA).

sequence and phylogeneTic analyses
The sequences were initially aligned using the default pair-wise and multiple alignment parameters in CLUSTAL X (Larkin et al., 2007) and adjusted m a n u a l l y u s i n g M e s q u i t e 3 . 2 ( M a d d i s o n & Maddison, 2017). Gaps were positioned to minimize nucleotide mismatches and treated as missing data in phylogenetic analyses. Poorly aligned regions that may interfere with phylogenetic estimation were located using trimAl v.1.2 with the -automated1 option (Capella-Gutiérrez, Silla-Martínez & Gabaldón, 2009) and excluded from the analyses.
Congruence of the datasets was assessed using a hierarchical likelihood ratio test (hLTR) implemented in Concaterpillar ver. 1.7.2 (Leigh et al., 2008).
Phylogenetic trees were inferred using the maximum likelihood (ML) method implemented in RAxML ver. 8.2.4 (Stamatakis, 2014) and the Bayesian inference (BI) method implemented in MrBayes ver. 3.2.6 x64 MPI (Ronquist et al., 2012). Several groups of accessions had identical sequences for some markers analysed. Each of these groups was represented by a single terminal in phylogenetic analyses. In the resulting trees, these terminals were added forming a polytomy composed of zero-length branches.
ML analyses included 200 searches starting from distinct randomized maximum parsimony trees. Branch support (BS) was evaluated based on 1000 rapid bootstrap replicates. To assess whether the number of replicates was sufficient, a posteriori bootstopping analysis was carried out with the extended majorityrule consensus tree as a criterion of convergence. All analyses were performed with GTR + G substitution model, which was most efficiently implemented and optimized for RAxML and thus recommended by the author (Stamatakis, 2014). For plastid data and combined nuclear and plastid data, partition schemes were inferred using PartitionFinder ver. 2.1.1 (Lanfear et al., 2016) with branch lengths linked and BIC options.
Bayesian analyses comprised two independent runs, each with four Monte Carlo Markov chains that were run for 40 000 000 generations with a sampling frequency of 1000 generations and default priors for the parameters of nucleotide substitution model. The initial 25% of saved trees were discarded and the results were summarized on the 50% majority rule consensus tree. The effective sample size (ESS) for the estimated parameters and the convergence of the independent runs were checked using Tracer ver. 1.6.0 (Rambaut et al., 2014). The model of nucleotide substitution for independent analyses of nrDNA ITS data was selected with ModelGenerator ver. 0.85 (Keane et al., 2006). For the analyses of plastid data and combined nuclear and plastid data, partition schemes and substitution models were inferred using PartitionFinder comparing the models implemented in MrBayes.
Two methods were used to identify problematic terminals, i.e. those that introduce putative topological conflict or decrease branch support: (1) compat.py (Kauff & Lutzoni, 2002 and (2) the RogueNaRok algorithm (Aberer, Krompass & Stamatakis, 2013). The first compares support values for trees inferred from different partitions (e.g. nrDNA and plastid DNA) and identifies terminals that occur in different and well-supported clades. We therefore compared clades occurring in nrDNA and plastid DNA trees choosing the values of BS ≥ 75 and PP = 1 as thresholds for ML and Bayesian analyses, respectively. In contrast, the RogueNaRok algorithm identifies a set of rogue terminals that, if pruned from the bootstrap trees, results in a consensus tree with more resolved polytomies (i.e. containing additional bipartitions) or with increased support values. We applied this algorithm to identify rogue taxa in ML analyses of combined data. The terminals identified by this algorithm were pruned from the ML bootstrap trees, which were then summarized on the best tree. Bootstrap support values for major clades were compared to those obtained from the analyses of unpruned trees.
To check whether the terminals identified with both methods introduced incongruence between nrDNA and plastid DNA, we excluded them from the matrices and repeated the hierarchical likelihood ratio tests (hLTR) with Concaterpillar. The data for this study have been deposited in TreeBASE, study number 21677 (http:// purl.org/phylo/treebase/phylows/study/TB2:S21677).

sequence characTerisTics
A l i g n m e n t o f t h e s e q u e n c e s wa s r e l a t i v e l y unambiguous. In the matrices trimmed with trimAl, the ITS region was more variable for the ingroup than the concatenated plastid DNA regions, with 25.5% variable positions as opposed to 10.0% in the latter (Table 1). Seventy-two ITS sequences were nonunique, distributed in 17 groups comprising two to 16 terminals each. In the concatenated plastid DNA data, 38 accessions were non-unique and distributed in 16 groups of two to four terminals each. However, in the matrix of combined nrDNA and plastid DNA data, only six pairs of accessions had all sequences identical. Each group of identical sequences was represented by a single terminal in subsequent phylogenetic analyses.
Hierarchical clustering of the markers based on likelihood-ratio tests concatenated all three plastid DNA regions with P = 0.043 and corrected α = 0.025; however, it failed to concatenate plastid DNA and ITS with P < 1 × 10 -6 and corrected α = 0.017. For ITS data, ModelGenerator with the Bayesian information criterion (BIC) selected the SYM + G nucleotide substitution model. PartitionFinder partitioned the data into ITS and plastid DNA with the SYM + G and GTR + G + I substitution models, respectively.

phylogeneTic analyses
Trees obtained from ML and Bayesian analyses were similar in topology; therefore, here we only show ML trees, whereas Bayesian trees are presented electronically in the Additional Supporting Information. We identify major clades with letters and their subclades with numbers; clades A-I refer to those recognized by Kurzyna-Młynik et al. (2008) and J-L are recognized here for the first time.
Separate analyses of nrDNA ITS and plastid DNA matrices of these data resulted in trees with notably different topologies. In the ITS trees, major clades of Ferulinae formed a basal polytomy (Figs 1 and S1), whereas in the plastid DNA trees relationships among those clades were better resolved, but also contained polytomies (Figs 2 and S2). However, several clades that were apparent in the ITS trees did not occur in the plastid DNA trees. Although some clades were present in the trees inferred from both datasets, numerous terminals were placed in different positions in these trees (Fig. S3). Compat.py identified 41 terminals that had significantly different positions in nrDNA and plastid DNA trees (Figs 1 and 2). Notably, the ITS tree was more congruent with the traditional taxonomic treatment of the genus than the plastid tree. For instance, all samples of F. gummosa Boiss. have identical ITS sequences and occurred in clade H, whereas each of them has a different plastid haplotype Table 1. Characteristics of the datasets used in phylogenetic analyses. Apart from sequence length variation, all numbers concern matrices without ambiguously aligned positions removed using trimAl program. Numbers in parentheses refer to the ingroup (Ferula)  1. Maximum likelihood tree inferred from analyses of nrDNA ITS sequence data for Ferulinae. Groups of accessions with identical sequences were represented in the analyses by single terminals and then resolved to polytomies with that was more closely related to representatives of other species than to conspecific individuals. Similar significant discrepancies occurred in F. ovina (Boiss.) Boiss. (clade A). Also noteworthy are the relationships among the four representatives of F. assa-foetida L. (clade H), although they were not supported by high BS and PP values and, therefore, these samples were not identified as conflicting by compat.py. Based on ITS data, terminals 0433 and 2051 are identical and were grouped with 0148, whereas sample 0359 was placed on a different branch (Fig. 1). The latter, however, had the same plastid DNA haplotype as sample 0148, whereas accessions 0433 and 2051 were placed on distant branches (Fig. 2). This example suggests that the discordance between nrDNA and plastid DNA data does not concern exclusively the terminals identified by compat.py as conflicting. Indeed, despite removing these 41 samples from the nrDNA and plastid DNA matrices, these datasets remained incongruent in Concaterpillar analyses (P < 1 × 10 -6 , corrected α = 0.017).
Major clades in the ML and Bayesian trees of combined data (Figs 3 and S4, respectively) were almost the same as those in the nrDNA trees, whereas the relationships among these clades were more similar to those in the plastid DNA trees. Leutea and Ferula were found to be sister taxa in Ferulinae; the monophyly of each genus was strongly supported both in ML and Bayesian analyses (BS ≥ 93%, PP = 1.0). In combined ML analyses of reduced data without conflicting terminals, major clades and their bootstrap support were similar to these in the tree obtained from full dataset (compare Figs 3 and S5).
In Ferula, there were 11 major clades in the combined ML tree (Fig. 3). All but clades H and L also occurred in nrDNA analyses (Fig. 1), whereas only four clades (B, D, L and J) were apparent in ML plastid trees (Fig. 2). In the Bayesian consensus tree from the combined data some clades collapsed into polytomies (Fig. S4).
The rogue taxa analysis using RogueNaRok identified 13 unstable terminals: six accessions belonging to clade A, one to clade E, one to clade G, four to clade H and F. kuhistanica Korovin, which does not have a definite placement (Fig. 3). The removal of these terminals from the dataset notably increased bootstrap support for clades E (from 69 to 87%) and H (from 12 to 49%); nevertheless, the bootstrap support for clade H remained low. The support for other clades changed only marginally (Fig. S5). Despite deleting these 13 taxa from the matrix, plastid DNA and nrDNA datasets remained incongruent in Concaterpillar analysis (P < 1 × 10 -6 , corrected α = 0.017).

DISCUSSION
incongruence beTween nrdna and plasTid dna daTa Several factors may have contributed to the detected incongruence between nrDNA ITS and plastid DNA data. We cannot conclusively exclude accidental contamination or human error at various stages of laboratory procedures. Due to a limited availability of plant material, samples were mostly collected from old herbarium specimens, sometimes from crumbled leaves or flowers stored in envelopes. Such materials may be accidentally contaminated with tissues of other species, e.g. those stored in the same herbarium cabinets. Additionally, some PCR primers may be more specific to the contaminant than to the nominal species. As we explained in the Material and Methods section, for some difficult DNA samples, the selected markers were amplified and sequenced in parts using various combinations of external and internal primers. Partial sequences were then assembled based on overlapping regions. These regions were usually evolutionarily conserved, particularly in the plastid DNA markers, and did not vary much among species. Therefore, a chimaeric contig from sequences originating from different species might have passed unnoticed. For the same reasons, a human error would not have been easily detected. When rechecking the data from our previous paper on Ferulinae (Panahi et al., 2015), we found indeed that one sequence of the rps16 intron (KJ660433) presumed to represent Leutea glaucopruinosa (Rech.f.) Akhani & Salimian was from a species of Orlaya Hoffm. (other plastid DNA markers from L. glaucopruinosa were correct and support the placement of this species in Leutea). Additionally, repeated sequencing for two species, F. tadshikorum Pimenov (rps16 intron KJ660415) and F. tuberifera Korovin (rpoB-trnC spacer KJ660729), indicated that the DNA samples were probably contaminated with other species; therefore, these accessions were not used in the analyses presented here.
Despite the aforementioned problems, there are good reasons to conclude that the discrepancies between ITS and plastid DNA data are not artefacts but have resulted from biological processes. First, we repeated PCR and sequencing for some samples that exhibited zero-length branches. Bootstrap support and posterior probability for nodes that also occurred in Bayesian 50% majorityrule consensus tree (see Fig. S1) are given along the branches. Accessions identified by compat.py as conflicting are marked with different colour for each group. Major clades are bracketed; clades A-I follow Kurzyna-Młynik et al. (2008). Outgroup taxa are omitted for simplicity.  striking incongruence and got the same results. Particularly, we confirmed that the three samples of F. gummosa with identical nrDNA ITS sequences have different plastid haplotypes that are more closely related to those occurring in other species than to each other. Second, the discordance between plastid DNA and nrDNA was observed for DNA samples processed in two laboratories: in Warsaw and in Tehran (the latter annotated 'K' in Table S1). Third, the discrepancies occurred more often in major clades than among them suggesting that they resulted from hybridisation or introgression between closely related species as exemplified by F. gummosa and its close relatives.
The evolutionary processes that may result in incongruence among molecular markers are hybridization, introgression, incomplete lineage sorting and homoplastic substitutions, particularly those that have a stabilizing effect on the stem-loop secondary structure in plastid DNA sequences. This incongruence may also result from sampling error of characters or/and taxa (Salichos, Stamatakis & Rokas, 2014). However, the recorded discrepancy is more likely biological than artefactual as exemplified by F. gummosa. Possible hybridization and introgression resulting in discrepancy between plastid DNA and nrDNA ITS data have also been reported in several other studies of Apiaceae (Lee & Downie, 2006;Spalik, Downie & Watson, 2009;Zhou et al., 2009;Bone et al., 2011;Yi, Jin & Wen, 2015;Banasiak et al., 2016).
Our results indicate that most intense reticulate evolution occurs among the species constituting clades A, E, G, H, I, K and L in ITS and combined data trees. These species along with members of clades B and C2 form a weakly supported superclade in plastid DNA trees (BS < 50%, PP = 0.99) comprising mostly Irano-Turanian species. Relationships in this superclade in the plastid DNA trees are fundamentally different from those occurring in the nrDNA trees. Conspecific samples more often group in the nrDNA trees than in the plastid DNA trees. The examples include F. assa-foetida, F. szowitsiana DC., F. gummosa, F. samarkandica Korovin, F. hedgeana Pimenov & Kljuykov etc. Moreover, the number of unique ITS sequences is rather low and several groups of species each have the same ribotype, particularly in clades A, G and H. This suggests that the exceptionally high number of species recognized in the genus is an artefact resulting from hybridization and taxonomic splitting rather than from diversification in the region. Intense reticulate evolution in Ferula explains why proposing an unambiguous hierarchical classification system of the genus is almost impossible. Hybridization and introgression blur boundaries among species and confound phylogenetic reconstructions, particularly when distant species exchange genetic material.
Newly described species may represent local hybrid swarms, and inadequate sampling, particularly from Central and West Asia, does not allow for a precise delimitation of taxa.

Towards a new classificaTion sysTem
There is no updated and comprehensive classification system of the genus. The most important regional and worldwide treatments of Ferula comprising infrageneric taxa include: (1) a revision for the Flora Orientalis (Boissier, 1872); (2) a synopsis of the genus in Engler and Prantl's Die natürlichen Pflanzenfamilien (Drude, 1898); (3) the taxonomic monograph of the genus (Korovin, 1947); (4) a regional treatment for Kazakhstan (Safina & Pimenov, 1984). Korovin's classification system of Ferula was adopted with some amendments in subsequent regional revisions for the former Soviet Union (Korovin, 1951), Turkey and the East Aegean islands (Peşmen, 1972) and the Flora Iranica region (Chamberlain & Rechinger, 1987). None of these treatments agrees with molecular phylogenetic trees (Kurzyna-Młynik et al., 2008;Panahi et al., 2015; this study).
In accordance with our inference of a highly reticulated scenario of evolution in Ferula, when proposing a new classification system of the genus, a strict concept of monophyly, particularly in the Irano-Turanian superclade, would be misleading. This superclade is monophyletic in plastid DNA trees and collapses into a polytomy in the ITS and combined trees (Fig. 3). We recognize it as subgenus Narthex (Falc.) Drude with eight sections corresponding to clades A, B, E, G, H, I, K and L. In contrast, clades J, D and C1 are supported by both datasets and constitute earlydiverging branches in plastid DNA trees. Accordingly, we recognize them as subgenera Sinoferula, Safinia and Ferula. In the latter, we also include clade C2, the affinity of which to C1 is supported by nrDNA and morphology. Below, we discuss these clades in the order from the bottom to the top of the combined tree (Fig. 3).

Clade J: subgenus Sinoferula
The monophyly of this small group is supported in all analyses (BS = 100% and PP = 1.0 in combined trees), but it has not hitherto been recognized as a natural taxon. Its included F. licentiana Hand.-Mazz., F. kingdon-wardii H.Wolff and F. olivacea (Diels) H.Wolff ex Hand.-Mazz. that are endemic to China and have one distinctive common feature: they are glabrous throughout, in contrast to all other Chinese congeners that are hispid, pubescent or puberulent at least on the abaxial side of leaves (She & Watson, 2005). Numerous Irano-Turanian congeners also lack indumentum (Chamberlain & Rechinger, 1987) but the geographical ranges of these groups do not overlap. The fruits of F. kingdon-wardii and F. olivacea have prominent dorsal and lateral ribs and single vallecular and two commissural vittae, which are septate (Wang et al., 2016). Similar fruits occur in F. bungeana Kigatawa, which was not included in our molecular analyses due to incomplete availability of plastid markers. Based on the nrDNA ITS sequence, this species is allied with clade A. The plant is, however, densely pubescent. Septate vittae also occur in clade C, mostly comprising Mediterranean taxa (Safina & Pimenov, 1990).

Clade D: subgenus Safinia
Although three species of clade D were identified by compat.py as introducing topological conflict, these discrepancies concerned within-group relationships and the clade was strongly supported in trees inferred from combined data (BS = 93%, PP = 1.0). It encompasses Central Asian endemics with the notable exception of F. hezarlalehzarica Ajani from Iran. Clade D has never been recognized as a separate taxon; however, it is reasonably well corroborated by fruit anatomy. Fruits of F. equisetacea Koso-Pol., F. koso-poljanskyi Korovin, F. grigorievii B.Fedtsch. and F. decurrens Korovin are characterized by a subcircular disposition of vittae, a development of a subepidermal layer of the ducts (these may, however, be obsolete in mature fruits of F. grigorievii), and the parenchymatic ribs, hypendocarp and funicle (Safina & Pimenov, 1990). However, this combination of characters is not unique in the genus as F. lipskyi Korovin is similar to these species (Safina & Pimenov, 1990;Safina, Ostroumova & Pimenov, 2014). Ferula lipskyi was placed in clade A based on ITS sequence only (Kurzyna-Młynik et al., 2008). Additional sampling from this species is necessary to confirm its distant position from the species of clade D. The similarity of F. hindukushensis Kitam. to F. koso-poljanskyi was noticed by Chamberlain & Rechinger (1987). When describing F. hezarlalehzarica from the Hezar and Lalehzar mountains in Iran, Ajani & Ajani (2008) indicated the two former species as its closest relatives. The name of this subgenus honours botanist Lucia K. Safina for her contribution to the taxonomy of the genus.

Clade C: subgenus Ferula
This group includes two distinct subclades, C1 and C2, that have each strong support from posterior probability and bootstrap analyses, but their sister relationship is supported only in nrDNA analyses. Seven species that form clade C1 are mostly Mediterranean, including F. communis L., the type species of the genus. This group has long been regarded as natural, although some distantly related taxa were also usually included. Safina & Pimenov (1990) investigated fruit anatomy of 16 representatives of subgenus Ferula sensu Korovin (1947) and concluded that F. communis, F. tingitana L., F. glauca L. and F. linkii Webb form a distinct group characterized by septate vittae, narrow and tangentially stretched elements of hypendocarp and a moderate sclerification of the pericarp. The most widely distributed member of this clade, F. communis, is diversified into local taxa that are often recognized as distinct species (Dettori et al., 2016).
Based on plastid data, F. stenocarpa, the only representative of clade C in Iran, is closer to other Irano-Turanian taxa than to its Mediterranean relatives. Therefore, it may be of hybrid origin. This species is morphologically similar to its congeners from clade C, and this affinity has been generally recognized (Boissier, 1872;Korovin, 1947). It differs from other members of this clade in having distinctly smaller fruits, 4 mm long, as opposed to 8-18 mm (Korovin, 1947;Chamberlain & Rechinger, 1987). Based on phylogenetic analyses of ITS data, two Anatolian endemics not included in this study, F. coskunii H.Duman & Sağıroğlu and F. mervynii Sağıroğlu & H.Duman, also belong to clade C (Elibol et al., 2012). Our preliminary analyses using ITS data place them close to F. stenocarpa (unpubl. data). These species have also smaller fruits than the Mediterranean congeners: 4-5 mm long in F. mervynii (Sağiroğlu & Duman, 2007) and 5-9 mm in F. coskunii (Duman & Sağiroğlu, 2005). We therefore recognize clades C1 and C2 as sections Ferula and Stenocarpa, comprising Mediterranean and Irano-Turanian taxa, respectively.

Clade L: section Glaucoselinum (Schischk.) Pimenov
The similarity between the two species in this clade was first noticed by Korovin (1962), who segregated them from Peucedanum L. into Talassia Korovin. Pimenov (1983) transferred them to Ferula and placed them in a separate section but without subgeneric assignment and repeated this treatment in the subsequent revision of the genus for Kazakhstan (Safina & Pimenov, 1984). Because the validity of Korovin's Talassia was problematic, Pimenov (1983) used a sectional name Glaucoselinum that was first introduced in Peucedanum.

Clade E: section Macrorrhiza Korovin
This group includes five species occurring in the Irano-Turanian region; all species are glabrous (Chamberlain & Rechinger, 1987). Its support in the combined ML analyses was rather weak (BS = 69%), which is not surprising because two inclusive species were placed in significantly distant positions in plastid DNA trees. Three representatives of this clade (F. oopoda (Boiss. & Buhse) Boiss., F. clematidifolia Koso-Pol. and F. varia Trautv.) have been subject to carpological studies and have similar mericarp anatomy identified as type VIII (Safina & Pimenov, 1984). However, superficially similar fruits occur in F. korshinskyi Korovin, a member of clade A.

Clade B: section Soranthus
Ferula karelinii Bunge and F. sibirica Willd. form a strongly supported clade in combined analyses (BS = 100%, PP = 1.0). These species were excluded from the genus by Korovin (1947) and recognized in monospecific genera Schumannia Kuntze and Soranthus Ledeb., particularly in the regional treatments for the former Soviet Union and Asia, even in some relatively recent accounts (Schischkin, 1951;Vinogradova, 2004;She et al., 2005). Their fruit morphology and anatomy lies, however, in the diversity range of Ferula (Pimenov & Kirillina, 1980). Both species have fruit type III based on the typology of Safina & Pimenov (1984); these authors placed F. karelinii in section Merwia (B. Fedtsch.) Koso-Pol., comprising species with fruit type II, while retaining F. sibirica in the monotypic section Soranthus.

Clade A: section Peucedanoides
This large clade encompasses species occurring in Central Asia and in the south-western part of the Irano-Turanian region. It includes all species of former Dorema examined for molecular markers (Panahi et al., 2015) and nearly all examined species of subgenera Peucedanoides and Dorematoides that are, however, intermingled. Its monophyly was reasonably well supported in combined analyses (BS = 88%, PP = 1.0), but most internal nodes received BS < 50%. Its two major subclades in the combined tree recieved no support in ITS and plastid DNA analyses.
Species formerly classified in Dorema have simple umbels that are arranged laterally on a flowering stem, whereas most umbellifers have compound umbels that terminate the growth of the stem. The umbels of species classified by Korovin (1947) in Dorematoides are proliferating: a new umbel grows in place of the central umbellule forming a pseudoverticillate inflorescence and restoring its monopodial growth. This type of inflorescence can be seen as transitory between terminal compound umbels characteristic for most Ferula spp. and lateral simple umbels occurring in former Dorema. Korovin (1939) attributed the name Dorematoides to Regel and Schmalhausen. In fact, this name was introduced by Regel (1878) as group 'E. Doremoides' in a key to some Ferula spp. (with letters identifying successive groups) and repeated in a similar context in a subsequent paper (Regel & Schmalhausen, 1878). 'Doremoides' included only F. schtschurowskiana Regel & Schmalh. ex Regel ('F. tschzurowskiana'). The name was unranked and without any authorship in contrast to some formal names listed in these papers. Moreover, in Regel (1878), a parallel group in 'Sectio II. Asa foetida' was named 'D. Juga vittata', subsequently changed to 'D. Jugivittatae' in Regel & Schmalhausen (1878), suggesting that these groups were intended as informal and hence not validly published. Korovin (1947) corrected the name to Dorematoides taking into account the declension of Greek noun dorema, dorematos ('gift') and provided a Latin diagnosis. Therefore, he should be regarded as the sole author of the name at the subgeneric rank. Such an interpretation solves the problem of its typification. Vinogradova (2004) lectotypified 'section Dorematoides Regel & Schmalh.' with F. caspica M.Bieb. However, neither Regel (1878) nor Regel & Schmalhausen (1878) included F. caspica in group 'Doremoides', but indicated F. schtschurowskiana as its sole element. Korovin (1947) placed six species in subgenus Dorematoides indicating the proliferating umbels as their major feature, which makes these species somewhat similar to Dorema. In molecular analyses, F. schtschurowskiana was placed in clade G, whereas F. caspica with two other members of subgenus Dorematoides [F. feruloides (Steud.) Korovin and F. dubjanskyi Korovin] were assigned to clade A that also includes former Dorema spp. The choice of F. caspica over F. schtschurowskiana better preserves the original meaning of Dorematoides sensu Korovin. Since Dorematoides was described by Korovin as a subgenus, the author of the combination at sectional rank is Vinogradova (2004), who published it inadvertently when citing Korovin's original name with full reference.

Clade K: section Pachycarpa
The core of clade K, F. gigantea B.Fedtsch. and F. trachyphylla Rech.f. & Riedl, was supported in all analyses. In ITS and combined ML analyses, F. diversivittata Regel & Schmalh. joined this group. Ferula gigantea and F. diversivittata were placed by Korovin (1947) in an unranked group (grex) Pachycarpa. Ferula trachyphylla was described from Afghanistan and regarded as closely allied to F. gigantea (Chamberlain & Rechinger, 1987), and this is confirmed by molecular data. We designate F. gigantea as the lectotype of Pachycarpa and use this name as the basionym for the sectional combination.

Clade I: section Euryangium (Kauffm.) Pimenov
The support for the monophyly of this group comes exclusively from nrDNA ITS data. Having been transferred to Ferula from Euryangium Kauffm., F. sumbul (Kauffm.) Hook.f. has usually been recognized in a separate subgenus (Drude, 1898), section (Pimenov, 1979) or informal grex (Korovin, 1947) with one or two presumed close relatives, the affinity of which has not been confirmed by molecular data. The fruit anatomy of F. sumbul is distinct (Safina & Pimenov, 1983); however, the other species from this clade have not yet been examined for fruit anatomical characters.

Clade G: section Scorodosma (Bunge) Boiss
Molecular data place 18 Central Asian and Iranian species in this clade. Although its monophyly is strongly supported in ITS and combined analyses (BS = 91%, PP = 1.0), the relationships among included species are poorly resolved. Many terminals were identified by compat.py as introducing topological conflict, and in the plastid DNA trees they were intermingled with some members of clade H. Conspecific samples of F. hedgeana and F. samarkandica grouped only in ITS analyses, whereas they were intermingled with other species in plastid and combined analyses. This clade includes F. foetida (Bunge) Regel, which had been regarded by Korovin (1947) and Chamberlain & Rechinger (1987) as having an isolated position in the genus and constituting a monospecific subgenus Scorodosma (Bunge) Drude. However, Safina & Pimenov (1984) suggested that its relatives include F. iliensis Krasn. ex Korovin and F. teterrima Kar. & Kir., the affinity of which is confirmed by our analyses.

Clade H: section Merwia
Species placed in clade H occur in the south-western part of the Irano-Turanian floristic region. Our results suggest that this group is subject to intense reticulate evolution blurring the boundaries among species. This is particularly troublesome because this group comprises several economically important species including F. assa-foetida and F. gummosa.
The identity of F. assa-foetida is contentious. A possible source of confusion is that the exudate asafoetida is harvested from several Ferula spp., not only representatives of clade H, e.g. F. alliacea Boiss., F. persica Willd. and F. narthex, but also distant relatives including F. foetida (Drude, 1898;Korovin, 1959). Linnaean F. assa-foetida was based on material collected by Kaempfer. Chamberlain & Rechinger (1987) compared the plate in Kaempfer (1712) and the type specimen in the Sloane herbarium at BM with the available herbarium material and identified only two gatherings that are conspecific with the type. One of these is Davis & Bokhari D.56275; this specimen was also used in our molecular studies (Number 0359). In nrDNA trees, it was placed close to F. pseudalliacea Rech.f., in agreement with Chamberlain & Rechinger (1987), who regarded F. assa-foetida and F. pseudalliacea as closely allied and possibly conspecific. Three other samples of F. assa-foetida grouped with F. alliacea, F. gabrielii and other species harvested for asafetida (Chamberlain & Rechinger, 1987).
The examined molecular markers provide conflicting results concerning the affinities of F. gummosa, a source of valuable galbanum. This oleo-gum-resin is subject to adulteration as the demand is higher than its production in Iran (Betti et al., 2010). Presumed adulterants include F. assa-foetida and F. ammoniacum (D.Don) Spalik et al. (Thomsen et al., 2004). Therefore, a clear circumscription of F. gummosa is also of considerable economic importance. In our analyses, all samples of F. gummosa have the same ITS sequence, which is also identical to those of F. badrakema Koso-Pol., F. linczevskii Korovin, F. undulata Pimenov & J.V.Baranova, and F. myrioloba Rech.f. In the plastid DNA tree, however, the representatives of F. gummosa were placed in three distant clades, sometimes with their companions from the ITS tree. These results suggest that interspecific boundaries in this group are unclear; therefore, presumed adulterants as assessed through phytochemical studies may reflect a hybrid origin of harvested populations. Chamberlain & Rechinger (1987) suggested a close proximity of F. gummosa, F. myrioloba and F. badrakema, which is confirmed by molecular trees. Pimenov & Kljuykov (1996) synonymized F. badrakema under F. gummosa and also reinstated F. galbaniflua Boiss. & Buhse in this group. Chamberlain & Rechinger (1987) regarded F. latisecta Rech.f. & Aellen and F. undulata as closely related if not conspecific, also in agreement with molecular data. However, they placed F. latisecta and F. gummosa in separate infrageneric divisions.

SUPPORTING INFORMATION
Additional Supporting Information may be found in the online version of this article at the publisher's web-site: Figure S1. Bayesian majority-rule consensus tree inferred from analyses of nrDNA ITS sequence data of Ferulinae. Groups of accessions with identical sequences were represented in the analyses by single terminals and then resolved to polytomies with zero-length branches. Posterior probabilities are given along the branches. Accessions identified by compat.py as conflicting are marked with a different colour for each group. Major clades are bracketed; clades A-I follow Kurzyna-Młynik et al. (2008). Outgroup taxa are omitted for simplicity. Figure S2. Bayesian majority-rule consensus tree inferred from analyses of combined plastid rps16 intron, rpoC1 intron and rpoB-trnC intergenic spacer sequence data of Ferulinae. See Figure S1 for details. Figure S3. Comparison of topologies of maximum likelihood trees inferred from nrDNA ITS data (left) and concatenated plastid DNA data (right). See Figures 1 and 2 for details. Figure S4. Bayesian majority-rule consensus tree inferred from analyses of combined nrDNA ITS and plastid rps16 intron, rpoC1 intron and rpoB-trnC intergenic spacer sequence data for Ferulinae. See Figure S1 for details. Figure S5. Maximum likelihood tree inferred from analyses of combined nrDNA ITS and plastid rps16 intron, rpoC1 intron and rpoB-trnC intergenic spacer sequence data for Ferulinae with the exclusion of 41 accessions indicated by compat.py as introducing topological conflict. Groups of accessions with identical sequences were represented in the analyses by single terminals and then resolved to polytomy with zero-length branches. Bootstrap support and posterior probability for nodes that also occurred in Bayesian majority-rule consensus tree are given along the branches. Table S1. Representatives of Scandiceae subtribe Ferulinae and outgroups examined in this study with respective GenBank reference numbers. Collector's data in bold face indicate the type specimens of respective species. Newly obtained sequences are identified with an asterisk behind the reference number.