Morphological and phylogenetic (nrITS barcode) analyses are conducted to clarify the taxonomic status of 10 Delphinium species (D. cuneatum, D. dyctiocarpum, D. duhmbergii, D. elatum, D. litwinowii, D. pubiflorum, D. puniceum, D. sergii, D. subcuneatum, D. uralense) grown 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 PCoA 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 independence of D. puniceum belonging to the Diedropetala section and D. elatum, D. uralense, D. dyctiocarpum and D. pubiflorum belonging to the Delphinastrum section. Finally, recursive partitioning is performed to develop the dichotomous key which can be used to differentiate between Delphinium species in the territory under study. In conclusion, we stress that the definition of species belonging to the Delphinastrum section is hindered by the presence of numerous intermediate or hybrid forms, on the one hand, and the impact of climatic conditions on the display of morphological (and specifically, taxonomically significant) traits, on the other.

Morphometric analysis

Morphometric analysis was conducted for 790 specimens of 46 populations of 10 Delphinium species (D. cuneatum, D. dyctiocarpum, D. duhmbergii, D. elatum, D. litwinowii, D. pubiflorum, D. puniceum, D. sergii, D. subcuneatum, D. uralense) grown in the south-east of European Russia (Astrakhan, Volgograd, Voronezh, Orenburg, Penza, Rostov, Samara, Saratov, Tambov and Ulyanovsk Provinces, and the Republics of Kalmykia, Mordovia and Bashkortostan) (Fig. 1, Table 2). Voucher specimens from each population were deposited in the SARBG herbarium (Botanical Garden of Saratov State University, Saratov, Russia). The relatively small number of populations subject to the research was due to their scarcity throughout the territory in question, especially in the steppe regions.

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, 2011). The generative growth stage was identified following Fyodorov’s recommendations (2003).

The analysis included from 10 to 30 specimens of each population. Some populations were represented by a smaller number of specimens due to the general small number of a population or the 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 display (ESM 1). The values of the quantitative parameters of a specimen’s flower, bracts, bracteoles and leaves were measured and listed in ESM 2 and 3.

The analyzed qualitative parameters considered significant for the differentiation between the species of the Delphinium genus (Tzvelev, 2001) were as follows: pubescence of stem, inflorescence axis, pedicles, sepals (inside and outside), ovary and fruits; color of nectar petals and, for the Diedropetala section, sepals; form 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) (Kashin et al., 2019).

For all statistic calculations, R environment was used. The whole dataset comprised 790 samples. The Shapiro-Wilk test was conducted to check the normality of the data. To analyze the dataset, several multidimensional approaches were applied. Non-metric multidimensional scaling (nMDS) (Ripley, 1996) with the Gower distance was carried out with all quantitative and qualitative parameters.

First, we performed the principal coordinate analysis (PCoA); 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.

Second, the Random forest analysis (Liaw, 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, Ripley, 2002) – another training method which analyzes the available classification to develop binary trees for the remaining dataset. The structure of such trees resembles a dichotomous key (Therneau et al., 2014). Therefore, we used the results of our recursive partitioning to generate the dichotomous key that can be used to define Delphinium species in the territory under study.

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.

10 non-coding and potentially highly-variable cpDNA regions: ndhC−trnV, ndhF, petA−psbJ, psbE−petL, trnT−psbD (Dong et al., 2012), trnC−petN (Lee, Wen, 2004) and trnL–trnF, trnH−psbA, trnS−trnG, matK (Jabbour et al., 2012) – were amplified with specimens of 10 Delphinium species from pure populations (containing no intermediate forms). The whole dataset was analyzed using the inter-gene transcribed spacer of ribosomal DNA (ITS1, 5.8S and ITS2) which was amplified with the NNC-18S10 and C26A primers (Wen, Zimmer, 1996).

The polymerase chain reaction (PCR) was carried out in 50 µl reaction volumes. The reaction mixture contained 10 µl ready-to-use PCR mix MaGMix (200 µM of each dNTP, 1.5 mM MgCl2, 1.5 U SmarTaqDNA-polymerase and buffer; Dialat Ltd., Moscow, Russian Federation), 35 µl of deionized water, 3.4 pmol of each primer and 5 µl of template DNA. The PCR was conducted in the Mastercycler gradient thermocycler (Eppendorf, Germany). The program was as follows: preliminary denaturation during 5 minutes at 95° C, then 35 cycles of 30 seconds each at 95° C, 30 seconds at 55° C and 2 minutes at 72° C, then final elongation during 10 minutes at 72° C. PCR products were cleaned with agarose gel and eluted with the NucleoSpin® Gel and PCR Clean-up kit (MACHEREY-NAGEL, Germany).

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 SINTOL 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 GenBank (accession numbers MT137556 - MT137651; MT154543-MT154615). The analysis also included sequences of several Delphinium specimens imported from the GenBank as outgroups.

The matrices were analyzed by the Maximum Likelihood method (ML) using the Tamura-Nei model in the MEGA ver. X program (Kumar, 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, 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 250000 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, 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 intra-group 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.

Some populations were represented by a smaller number of specimens due to the general small number of a population or the 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%.