Taxonomic revision of Delphinium (Ranunculaceae) in the South-East of European Russia

Morphological and phylogenetic (nrITS barcode) analyses are conducted to clarify the taxonomic status of 10 Delphinium species (D. cuneatum, D. dictyocarpum, D. duhmbergii, D. elatum, D. litwinowii, D. pubiflorum, D. puniceum, D. sergii, D. subcuneatum, D. uralense) occurring in the South-East of European Russia. The morphometric analysis is carried out with 22 quantitative and 32 qualitative parameters. Based on all parameters, the nMDS supports the differentiation of D. puniceum, D. sergii, D. uralense and D. pubiflorum whereas other studied taxa remain undistinguished. Furthermore, the Random forest analysis and the MrBayes phylogenetic analysis of ITS sequences confirm the species status of D. puniceum belonging to D. sect. Diedropetala and D. elatum, D. uralense, D. dictyocarpum and D. pubiflorum belonging to D. sect. Delphinastrum. Finally, recursive partitioning is performed to develop the dichotomous key which can be used to identify the species of Delphinium in the study area. In conclusion, we stress that the recognition of species belonging to D. sect. Delphinastrum is hindered by the presence of numerous intermediate or hybrid forms, on the one hand, and the impact of climatic conditions on the morphological (and, specifically, taxonomically significant) traits, on the other hand.


Introduction
Delphinium L. is a genus of primarily perennial herbaceous plants in Ranunculaceae. The genus is comprised of over 300 species occurring in the temperate and subtropical zones of the Northern hemisphere as well as in the high mountain regions of Western and Eastern Sub-Saharan Africa (Tzvelev 2001;Jabbour and Renner 2012;Chartier et al. 2016).
In the past twenty years, a series of molecular studies have shed some light on the phylogeny of Ranunculaceae. However, relationships and limits between the family's subfamilies and tribes remain problematic (Cossard et al. 2016). At lower ranks, relationships between taxa are even more uncertain. For instance, Delphinieae Schroed. (Aconitum L., Delphinium L., Consolida (DC.) S.F.Gray and Aconitella Spach), one of the largest tribes in Ranunculaceae (Jabbour et al. 2009), is still rather understudied, which leads to controversies in understanding of the scope of the perennial Delphinium and the status of the annual Consolida and Aconitella. The available molecular 1 3 59 Page 2 of 15 phylogenetic studies are rather in favor of the synonymous status of the latter two genera Renner 2011, 2012).
Furthermore, there is considerable uncertainty regarding the taxonomic structure of Delphinium s.s., e.g., the phylogenetic relationship between D. subgen. Delphinium and the monotypic East Asian endemic D. sect. Anthriscifolium (Xiang et al. 2017). The problem is further complicated by the fact that new Delphinium species continue to be described from the Urals (Kulikov 2000), the Altai Mountains (Ebel 2007), Siberia (Ebel 2007;Kurbatskiy 2017) and the Caucasus (Luferov 2011), which need to be taxonomically assessed.
There have been several attempts to clarify the phylogenetic structure of taxa in Delphinium. First, Koontz et al. (2009)  Gray). However, the maximum parsimony analysis run by the authors also allowed them to conclude that within the consensus phylogenetic tree generated from DNA sequences of 62 out of 67 studied species, the majority of taxa comprising D. sect. Diedropetala have insignificant bootstrap support and small resolution.
Also, molecular markers were used to investigate the phylogenetic relations in D. sect. Delphinastrum Renner 2011, 2012). Renner (2011, 2012) examined specimens of 54 out of 150 section species including two species subject to our research (one herbarium specimen per each species), i.e. D. elatum L. (from the Urals) and D. cuneatum Steven ex DC. (from Samara Province). Upon the analysis of plastid and nuclear DNA sequences (ITS-1-5.8S-ITS-2 and trnL-trnF) by the maximum likelihood method, D. elatum and D. cuneatum were not differentiated, although they are believed to belong, according to Malyutin (1987), to two different subsections of D. sect. Delphinastrum or, according to Tzvelev (2001), to two different species aggregates. In Flora Europeae (Tutin et al. 1993), D. elatum and D. cuneatum are also treated as two separate species. Table 1 shows the interpretation of Delphinium taxa in Eastern Europe by various authors (Malyutin 1987;Tutin et al. 1993;Tzvelev 2001; https:// wcvp. scien ce. kew. org/). Thus, given the conflicting views on the phylogeny of the genus in question, the objective of the present research is to define the taxonomic status and phylogenetic structure of Delphinium in the South-East of European Russia with the help of morphological and molecular methods.

Morphometric analysis
As a working hypothesis, we accepted Tzvelev's understanding (2001) of the taxonomic structure of the genus.
In order to maintain homogeneity in morphological data assessment, we limited the analysis to mature plants (Sharma and Pandit 2011). The generative growth stage was identified following recommendations of Fedorov (2003).
The analysis included from 10 to 30 specimens from each population. Some populations were represented by a smaller number of specimens due to a small size of the population or a small proportion of mature specimens. Whenever a specimen had a gap in morphometric data in the final matrix, the missing quantitative value was filled in by the mean variable value. The proportion of missing data was below 1%.
The analysis involved 54-32 qualitative and 22 quantitative morphometric parameters. As for the qualitative parameters, values were assigned to each specimen based on a parameter code (Online Resource 1). The values of the quantitative parameters of a specimen's flower, bracts, bracteoles and leaves were measured and listed in Online Resources 2 and 3.
The analyzed qualitative parameters considered taxonomically diagnostic in Delphinium (Tzvelev 2001) were as follows: pubescence of stem, inflorescence axis, pedicles, sepals (inside and outside), ovary and fruits; color of nectar petals and, in D. sect. Diedropetala, sepals; shape of blade base; divergence of lowest blade lobes; gap between lowest blade lobes; degree of leaf dissection; shape of bracts and bracteoles. In addition to the listed qualitative parameters, we analyzed a number of other qualitative (pubescence of spur and petiole) and quantitative parameters (plant height; shoot length; number of shoot leaves; number of paraclades; length of blade; length of non-dissected part of blade; width of blade; width of base of middle segment of central lobe; width of base of central lobe) . For all statistic calculations, R environment (R Core Team 2020) was used. The whole dataset comprised 790 samples. The Shapiro-Wilk test was conducted to check the normality of the data. To analyse the dataset, several multidimensional approaches were applied. Non-metric multidimensional scaling (nMDS) with the Gower distance was carried out with all quantitative and qualitative parameters (Venables and Ripley 2002;Dragulescu and Arendt 2020;Shipunov et al. 2020). First, the non-metric multidimensional scaling (nMDS) with the Gower distance was carried out with all quantitative and qualitative parameters (Venables and Ripley 2002;Dragulescu and Arendt 2020;Shipunov et al. 2020). Hulls were constructed. The distances between the hull centers as well as the degree of hull overlap were measured. Also, morphometric and geographic data were subject to the correlation analysis with the nMDS axes (Schnute et al. 2019).
Second, the Random forest analysis (Liaw and Wiener 2002) was applied to all data subsets; the earlier predicted species identities were taken into account. For the analysis, we used the "classification and training" method described in Shipunov and Efimov (2015).
Third, and finally, we conducted recursive partitioning (Venables and Ripley 2002), i.e. another training method which analyses the available classification to develop binary trees for the remaining dataset. The structure of such trees resembles a dichotomous key (Therneau and Atkinson 2019). Therefore, we used the results of our recursive partitioning to generate the dichotomous key that can be used to identify Delphinium species in the study area.

Molecular analysis
DNA was extracted from the perianth leaves and leaflets previously dried in silica gel. DNA was extracted using the NucleoSpin® Plant II kit (MACHEREY-NAGEL, Germany) according to the manufacturer's protocol.
Sequencing was performed in the ABI PRISM 3130 XL sequencer using the BIG DYE TERMINATOR kit ver. 3.1 according to the manufacturer's protocol from the SIN-TOL Company (Moscow, Russian Federation).
Forward and reverse sequences were edited and manually aligned in the BioEdit 7.0.5.3 program (Hall 1999). The obtained DNA sequences were submitted to the Gen-Bank (Online Resource 4). The analysis also included sequences of several Delphinium specimens imported from the GenBank.
The matrices were analyzed by the Maximum Likelihood method (ML) using the Tamura-Nei model in the MEGA ver. X program (Kumar et al. 2018). The evolutionary direction of sequence changes was inferred by the outgroup comparison. All most parsimonious trees were summarized by the strict consensus method. Bootstrap analyses (1000 replicates) were performed to assess support of the clades (Felsenstein 1985). Bayesian phylogenetic analyses were performed with MrBayes 3.1.23 (Ronquist and Huelsenbeck 2003). The sequence evolution model was evaluated by the Akaike Information Criterion (AIC) in jModeltest 3.7 (Darriba et al. 2012). Two independent analyses with four Markov chains were run for 10 million generations, sampling a tree every 100 generations. 25% of initial trees were discarded as burn-in. The remaining 250,000 trees were combined into a single data set, and a majority-rule consensus tree obtained. Bayesian posterior probabilities were calculated for that tree in MrBayes 3.1.23.
Also, the statistical parsimony analysis was performed according to the algorithm described in Templeton et al. (1992) and carried out in TCS v. 1.21 (Clement et al. 2000). In general, the method evaluates the unrooted haplotype network, and all the links between haplotypes with  95% probability. Only nucleotide substitutions were taken into account; indels were treated as the missing data. The analysis of molecular variance (AMOVA) was conducted in the ARLEQUIN 3.5 program (Excoffier and Lischer 2010) based on the pairwise comparison of specimens with 1000 permutations. To evaluate the dispersion between the species and species groups, we used the F CT (inter-group dispersion) and F SC (inter-population intragroup dispersion) values which were to be maximum and minimum respectively. The correlation between genetic and geographical distances between populations was evaluated by the Mantel test.

nMDS of morphological characters
The nMDS analysis with 54 morphological traits of 790 specimens of 10 Delphinium species resulted in identifying the clouds of points representing specimens of each species (Fig. 2a). In the top left corner of the axes space, there are three detached hulls corresponding to the specimens of D. puniceum, D. sergii and D. uralense. The D. pubiflorum hull is rather separated as well. In the top right  Fig. 2b. 12 specimens that do not join any hulls fall into the "outliers" category (see Online Resource 5). The correlation analysis with the nMDS axes revealed that the following parameters have the highest loading in the first axis (V1): color of nectar petals; shape of bracts; pubescence of ovary and fruits; degree of leaf dissection; pubescence of stem; pubescence of pedicles; pubescence of leaf blade; shape of bracteoles; disposition of bracteoles; number of shoot leaves; width of base of central lobe; and width of base of middle segment of central lobe. In the second axis (V2), shape of leaf blade and inner pubescence of sepal have the highest loadings (Fig. 3a, Online Resource 5).   Fig. 3 Correlations between axes of non-metric multidimensional scaling and: a measurement characters, b geographic coordinates Finally, the nMDS analysis with geographical coordinates detected a positive correlation between the first axis (V1) and latitude and a negative correlation between the second axis (V2) and longitude (Fig. 3b).

Random forest and recursive trees
In the Random forest analysis, we used our set of "species" as a training dataset. Upon the training, the Random forest method was applied to the whole dataset comprising 790 specimens. As a result, the species identities were predicted based on the identity descriptions fed during the training (see Online Resource 5). The rate of classification error amounted approximately to 7%.
Based on all morphological (quantitative and qualitative) parameters and taking into account the predicted species identities, we built two binary trees (Fig. 4). The first tree was generated in accordance with Tzvelev's (2001) interpretation of the studied species (Fig. 4a), while the second tree united D. cuneatum, D. duhmbergii, D. litwinowii and D. subcuneatum in one species with a priority name D. cuneatum (Fig. 4b). Considering that, of the two trees, the second one drew a clearer and simpler distinction between the studied taxa, it was selected as a premise for the subsequent generation of the dichotomous key.

Molecular analysis
The sequencing of 10 cpDNA regions gave the results with extremely low resolution. Singular (most likely occasional) substitutions (1 − 2 substitutions per region) were found in 8 out of 10 examined regions (see Online Resource 4). Therefore, the subsequent analysis was limited to the ITS region.
The statistical parsimony analysis in the TCS program subdivided all specimens into 3 evolutionary networks (Online Resource 6). The first network is comprised of D. puniceum specimens, the second -of D. elatum and D. dictyocarpum specimens and the third -of all the remaining specimens. The third network has a web structure with numerous links between haplotypes.
On the phylogenetic tree built by the MrBayes method (Fig. 5), the specimens of 9 studied species (excluding D. sergii) are subdivided into two clusters. The first cluster includes all D. puniceum specimens. The second cluster is comprised of three first-order subclusters. The first subcluster consists of the specimens of D. dictyocarpum population no. 1 occurring on the border between Orenburg and Samara Provinces, one D. elatum specimen from the Urals imported from the GenBank (JN573517), the specimens from D. pubiflorum populations no. 3 and 7, and one specimen from D. pubiflorum population no. 1.
The second subcluster is subdivided into two secondorder subclusters: one comprising all specimens of the Northernmost D. dictyocarpum population and the other comprising the specimens of D. elatum populations located in the Urals.
The third subcluster consists of the following secondorder subclusters. The first second-order subcluster contains the specimens of D. dictyocarpum populations no. 3 and 6 -both occurring in the Southern Urals to the South of the above-mentioned D. dictyocarpum population. The second second-order subcluster is subdivided into the following third-order subclusters. The first contains the specimens of all D. uralense populations. The second is subdivided into two fourth-order subclusters: one comprising all the remaining D. pubiflorum specimens and the other comprising the specimens of all D. agg. cuneatum species (D. litwinowii, D. cuneatum, D. subcuneatum, D. duhmbergii) and the specimens of three Western D. elatum populations (from the East of Orenburg Province and the North of Samara Province) without any differentiation between the species.
The AMOVA shows that of all the proposed taxonomies, three are the most feasible (Table 3). In the first taxonomy, specimens are joined into taxa (following the conventional approach in the Russian botanical literature). In the second taxonomy, they are united into groups, or aggregates (following Tzvelev's (2001) interpretation). In the third taxonomy, specimens are grouped according to our second binary tree built by the Random forest. In the latter case, we observe the minimum dispersion between population and the maximum dispersion between groups.
Finally, the Mantel's test does not reveal the statistically reliable correlation between the F ST (the index of fixation between populations in a group) pairwise values matrix and the matrix of pairwise geographical distances between populations (r = 0.110944; p = 0.069000).

Discussion
The nMDS analysis based on all quantitative and qualitative parameters differentiate between D. puniceum, D. sergii, D. uralense and D. pubiflorum, but the other taxa remain undistinguished. In our view, this confirms the complexity of phylogenetic and biogeographical relationships between D. elatum, D. dictyocarpum, D. cuneatum, D. subcuneatum, D. duhmbergii and D. litwinowii. Also, the nMDS analysis strongly supports the differentiation between D. dictyocarpum, D. pubiflorum and D. uralense which, according to Tzvelev (2001), belong to D. agg. dictyocarpum. Figure 1 shows that, geographically, these taxa are also well separated from each other. Delphinium pubiflorum is confined to the West, D. dictyocarpum -to the North-East, and D. uralense -to the South-East of their group range with virtually no geographical overlap between them. Notably, of the three taxa, D. uralense is the most isolated in the nMDS axes space and, being the most xerophytic in D. agg. dictyocarpum, the closest to even more xerophytic D. puniceum and D. sergii. Nevertheless, D. puniceum and D. sergii are still rather remote from the above-mentioned taxa.
As for D. elatum, in the nMDS axes space, it is separated not from D. dictyocarpum but from the other D. agg. dictyocarpum taxa. Geographically, D. elatum also has a vast contact zone with D. dictyocarpum in the North-East and East of their group range, whereas there is no geographical contact between D. elatum and the other D. agg. dictyocarpum species.
As Fig. 2 shows, the species of D. agg. cuneatum are the most difficult to distinguish in the nMDS axes space. These species are not differentiated by the nMDS and Random forest analyses either. In the nMDS axes space, they form a joint point cloud which overlaps with the D. elatum and D. dictyocarpum hulls (Fig. 2).
Apparently, this counts against, first, the separate status of the taxa belonging to D. agg. cuneatum and, second, the  30 D. cuneatum, D. litwinowii, D. duhmbergii, D. subcuneatum vs. D. dictyocarpum, D. pubiflorum, D. uralense 18 D. cuneatum, D. duhmbergii, D. litwinowii, D. pubiflorum, D. subcuneatum vs. D. uralense, D. dictyocarpum, D. 89 D. cuneatum, D. dictyocarpum, D. duhmbergii, D. litwinowii, D. pubiflorum, D. subcuneatum, D. uralense D. elatum from D. dictyocarpum in D. agg. cuneatum (although, according to Tzvelev (2001), D. elatum and D. dictyocarpum belong to different species aggregates). However, if one turns to the geographical distribution of all aforementioned taxa, it becomes clear that the D. cuneatum range overlaps with the D. elatum range and is adjacent to the Western part of the D. dictyocarpum range. Likewise, the D. litwinowii range is adjacent to the D. pubiflorum and D. duhmbergii ranges. Therefore, as Fig. 2 shows, in the nMDS axes space, it is the D. litwinowii hull that forms a joint hull between the D. agg. cuneatum taxa, on the one hand, and closely approaches the D. pubiflorum hull, on the other. As for D. dictyocarpum and D. elatum, their joint hull is formed by the specimens of both species occurring in the West of the species contact zone.
Thus, we propose that complex relationships between the taxa in D. sect. Delphinastrum may be attributed to the geographic and ecotopic differences of these taxa brought about by the influence of climatic conditions. This hypothesis is confirmed by the results of our analysis of the correlation between the studied morphological parameters and the geographical coordinate gradients (Fig. 3b).
Furthermore, this hypothesis is supported by the results of our earlier study of 16 populations of 5 Delphinium taxa (D. puniceum, D. sergii, D. pubiflorum, D. duhmbergii and D. litwinowii) in the Lower Volga region and the adjacent territory . The numeric analysis with 61 qualitative and quantitative parameters revealed that, in the nMDS axes space, the distribution of clouds representing specimens of various Delphinium populations and species corresponds to the continentality gradient, while the degree of hull overlap greatly depends upon the weather conditions of a certain year.
In this regard, we propose that the obtained results, namely the high morphological variability of Delphinastrum taxa, the morphological overlap between them, the impact of climatic conditions on the display of morphological parameters in specimens and the presence of outliers in their nMDS axes space, are pointing not only to the probable hybrid nature of the studied specimens but also to the instability and genetic variability of many of their taxonomically significant traits. This might explain the notorious difficulty in species-level identification of Delphinium herbarium specimens as well as the frequent re-evaluation of their species identity .
In the Russian botanical literature, there have been several studies concerning the role of climatic conditions in the inter-species and inter-population variability of Delphinium taxa (Malyutin 1987(Malyutin , 1992Fedorov 2003). For instance, the degree of the ecotopic and geographic variation of Delphinium species were examined in Fedorov (2003) within D. agg. elatum, D. agg. dictyocarpum and D. agg. cuneatum occurring in the Southern Urals. In 20 population samples of D.
agg. elatum alone, Fedorov (2003) identified 81 variants of pubescence (of various types and degrees) representing phenotypes of 3 species and 5 varieties belonging to the studied aggregate. According to Fedorov (2003), under moderate climate conditions (plain broadleaf and mixed forests), tall non-pubescent D. elatum var. elatum plants with wide leaf blades were widespread. Under colder and more humid conditions, pubescent phenotypes and varieties of D. elatum (D. elatum var. glandulosum) and D. alpinium species occurred.
According to Fedorov (2003), the spatial distribution of D. agg. elatum species, varieties and phenotypes showed their confinement to particular altitudes and latitudes as well as to open or forested habitats. A number of samples collected on the borders between species or variety ranges exhibited a high morphological variability. However, in one sample from a mixed population comprising two D. elatum varieties (D. elatum var. elatum and D. elatum var. glandulosum), two D. alpinum varieties (D. alpinum var. alpinum and D. alpinum var. hebecarpum) and D. nurguschense species, the primary factor defining the morphological diversity was terrain.
In our molecular analysis of the ITS region, all the employed methods give a sufficient support for differentiating D. puniceum as a separate species. Other specimens are subdivided by the statistical parsimony method into two evolutionary networks (Online Resource 6): one network comprising all Uralian D. elatum and D. dictyocarpum specimens and another network comprising all the other specimens. The fact that the second network has the web structure with the links between haplotypes confirms their hybrid nature and the existence of numerous intermediate forms.
According to the MrBayes analysis, D. agg. cuneatum species are not separated from one another. The fact that this cluster contains D. elatum specimens from populations no. 1, 2 and 3 located on the Western edge of the species range (Samara and Orenburg Provinces) confirms that these populations have the intermediate or hybrid nature. The same goes for the specimens of one Western D. dictyocarpum (DC-1) population, the haplotype of which is close to the haplotype of the specimens of three Southern D. pubiflorum populations (P-1, P-3 and P-7). Thus, in the phylogenetic tree, both D. elatum and D. dictyocarpum form a well separated cluster as well as the haplotypes of intermediate or hybrid forms occurring on the Western margins of their ranges. As for D. dictyocarpum, the haplotypes of intermediate forms are similar to the haplotypes of Southern D. pubiflorum specimens. As for D. elatum, they are close to the haplotypes of specimens belonging to various D. agg. cuneatum species.
Furthermore, in the phylogenetic tree, the haplotypes of Northern D. pubiflorum populations are also close to the haplotypes of all D. agg. cuneatum species, populations of which are adjacent to D. pubiflorum ones. The specimens of D. dictyocarpum population no. 2 are close to the specimens of the adjacent D. elatum populations occurring in the Urals. The specimens of D. dictyocarpum populations no. 3 and 4 are close to D. uralense specimens; these populations are geographically the most proximate as well. Delphinium uralense and D. puniceum are well separated from all the other species.
Thus, the phylogenetic tree of the ITS sequence haplotypes built by the MrBayes method proves the separate status of D. puniceum. Taking into account the results of the morphological analysis, the separation of D. sergii is also beyond doubt. In D. sect. Delphinastrum, D. elatum, D. pubiflorum, D. dictyocarpum and D. uralense should be recognized as separate species as well. At the same time, D. elatum, D. pubiflorum and D. dictyocarpum have intermediate or hybrid specimens found on the Western margins of their species ranges. Finally, D. litwinowii, D. cuneatum, D. subcuneatum and D. duhmbergii should be considered synonyms with the priority name D. cuneatum.

Conclusions
The results of our study confirm that, in the South-East of European Russia, the majority of the studied Delphinium taxa have the species status. The exceptions are D. litwinowii, D. cuneatum, D. subcuneatum and D. duhmbergii which are synonyms with the priority name D. cuneatum.
Our results also point to the difficulties hindering the definition of the taxonomy of Delphinium in the study area: namely, the presence of numerous intermediate or hybrid forms as well as the impact of climatic conditions on the expression of morphological -and specifically, taxonomically significant -parameters in Delphinium.
We believe that to resolve these difficulties, we need, first, to expand the geographical range and taxonomical scope of research by including the species from the Northern Caucasus and Siberia and, second, to conduct the whole genome sequencing of the species in D. sect. Delphinastrum.
The following key can be used to differentiate between Delphinium species in the study area: