Three Tertiary Euphorbia species persisted in the forests of the Balkan Peninsula

The Balkan Peninsula is renowned as a refugium of forest species, which were more widespread in the Tertiary. We here investigated three species from one of the largest genera of flowering plants, Euphorbia, which survived the Pleistocene glaciations in the Balkan Peninsula, but responded differently to Holocene warming. Using ITS sequences, multivariate morphometrics and relative genome size (RGS) measurements, we explored relationships among E. amygdaloides, E. heldreichii and E. orjeni, for which different taxonomic treatments have been proposed. The ITS-inferred phylogenies indicate that all three species form independent evolutionary lineages, which have similar RGS, but differ morphologically, and should thus be treated as independent species. The most enigmatic of them, E. orjeni, was found for the first time after 100 years at its type locality in Montenegro, and, based on herbarium revision, was discovered also in the surroundings of Belgrade. Euphorbia heldreichii is more widespread, distributed in the southern Balkan Peninsula, whereas E. amygdaloides occurs in forests throughout Europe, North Africa and Western Asia. Our study emphasises the importance of further botanical exploration of the Balkan Peninsula, a hotspot of European biotic diversity.


Introduction
The current distribution of biota throughout the Northern Hemisphere has been greatly affected by the climatic fluctuations and changing geography of the last 65 million years, i.e. during the Tertiary and Quaternary periods (Wolfe 1975;Milne and Abbott 2002). Aridification and cooling over the last 15 million years drastically changed the vegetation of Europe, where species of previously widespread Tertiary humidtemperate forests found their refugia in the Balkan Peninsula, Carpathians and adjacent Anatolia (Milne and Abbott 2002; Lendvay et al. 2016). With the decline of evergreen laurel forests, deciduous forest communities as well as evergreen gymnosperm forests expanded, but their ranges shrunk again during Pleistocene glaciations (Tiffney 1985) and modern Tertiary relict floras comprise mostly deciduous woody taxa and a small number of herbs (Milne and Abbott 2002).
The Balkan Peninsula is renowned as a refugium of several, previously more widespread, tree and shrub species. Aesculus hippocastanum L., Forsythia europaea Degen & Bald. (for which it is not clear if previously had larger distribution) and Picea omorika (Pančić) Purk. are examples of species that nowadays have very restricted distributions (Xiang et al. 1998;Aleksić and Geburek 2014;Ha et al. 2018), whereas Pinus heldreichii Christ, Syringa vulgaris L. and Tilia tomentosa Moench (Li et al. 2012;Saladin et al. 2017) cover larger areas on the Balkan Peninsula and beyond. In addition, many other tree species that in the Holocene expanded their ranges north of the Alps, had their glacial refugia on the Balkan Peninsula (Petit et al. 2002;Grivet and Petit 2003;Magri et al. 2006), which enabled persistence and diversification of understory forest species in this part of Europe (Kramp et al. 2009;Willner et al. 2009;Slovák et al. 2012;Rešetnik et al. 2016).
Euphorbia L., one of the three largest genera of flowering plants, exhibits a large diversity of growth forms (Horn et al. 2012). Its Eurasian radiation, E. subgen. Esula Pers. including about 480 species, diversified mostly in open, forest-free areas, especially in warm and dry habitats of the Mediterranean and Irano-Turanian regions (Riina et al. 2013). One of the exceptions is E. sect. Patellares (Prokh.) Frajman, which includes several mesophilous forest perennials, characterised by raylet leaves connate at the base (Riina et al. 2013). This lineage originated in the early Miocene, 17.7 Ma, but started to diversify only at the boundary between Miocene and Pliocene 5.6 Ma (Horn et al. 2014), giving rise to 14 species (Riina et al. 2013).
The most widespread of them is the European forest understory species E. amygdaloides L., distributed from Great Britain and Spain in the west over most of Central and Southern Europe to North Africa (Algeria, Tunis) and Western Asia, where it extends its range to the Caspian Sea area in the east (Radcliffe-Smith and Tutin 1968;Meusel et al. 1978;Radcliffe-Smith 1982). In addition, E. heldreichii Orph. was described from the southern Balkan Peninsula, where it grows in deciduous and pine forests in Albania and Greece (Radcliffe-Smith and Tutin 1968;Aldén 1986;Barina 2017). Aldén (1986) argued that it has seeds and capsules identical to E. amygdaloides and differs from it "only in having all the ray leaves and rays in whorls" and should thus be treated as parapatric subspecies E. amygdaloides subsp. heldreichii (Orph.) Aldén, which was followed also by Govaerts et al. (2000) and Dimopoulos et al. (2013).
Another Balkan endemic species similar to E. amygdaloides was described by Beck Mannagetta (1920) from Mt. Orjen in Montenegro (Veselý 1890) and to our knowledge remained known only from the type specimen until now. The species was ignored by Radcliffe-Smith and Tutin (1968), but also by Rohlena (1942) and Pulević (2005), who provided the most recent comprehensive treatments of the flora of Montenegro. On the other hand, Govaerts et al. (2000) included it in synonymy with E. amygdaloides and Riina et al. (2013) recognised it as a valid species, but it remained unclear whether this taxon is just a richly branched form of E. amygdaloides, as could be hypothesised based on the type material (B. Frajman, personal observation), or it deserves independent taxonomic status. Moreover, it is unknown where this species is distributed.
Taking into account uncertainties regarding the status of and relationships among widespread E. amygdaloides and the Balkan endemics E. heldreichii and E. orjeni, our aim is to explore phylogenetic relationships among these taxa and other co-occurring members of E. sect. Patellares and to investigate the molecular and morphological differentiation as well as possible relative genome size differences among the three before-mentioned taxa. To achieve this, (1) we searched for the existence of E. orjeni at its type locality and explored whether it is present elsewhere in the Balkan Peninsula by revising different herbarium collections. (2) We further inferred phylogenetic relationships among the taxa using nuclear ribosomal ITS sequences and (3) performed a morphometric analysis of all three taxa to elucidate which morphological characters potentially differentiate among them. (4) Using flow cytometry, we investigated variation in relative genome size in the three taxa, as 2n = 20 has been established for E. amygdaloides and E. heldreichii, but for E. amygdaloides occasionally also 18 chromosomes have been reported (Rice et al. 2015 and references therein). Finally, (5) we provide a revised taxonomic treatment for all three taxa and provide an identification key.

Plant material
Plant material of E. amygdaloides, E. heldreichii and E. orjeni for RGS estimation, molecular and morphometric analyses was collected in the field between 2005 and 2019. Euphorbia orjeni was searched for at its classical locality by taking the path described by Veselý (1890), who collected the only known herbarium specimen of this taxon. Molecular and RGS analyses were based on silica gel-dried leaf material, whereas morphometric measurements were performed on herbarium vouchers. In addition, we used four ITS sequences (one accession of E. amygdaloides and three of outgroup taxa) from previous studies (Frajman and Schönswetter 2011;Riina et al. 2013) deposited in GenBank and sequenced eight additional ingroup herbarium vouchers (deposited at BEOU and W) as well as six outgroup accessions of E. characias L. and E. macroceras Fisch. & C.A.Mey. (Online Resource 1). Moreover, 20 additional herbarium specimens from the herbaria W and WU were used for morphometric analyses. In total, 56 populations of E. amygdaloides (41), E. heldreichii (12) and E. orjeni (3) were studied, the numbers reflecting the differences in species' distributions. In most cases, the populations were represented by a single specimen in ITS and morphometric analyses, as we were not interested in intra-population variability, but in differences among the three taxa. In RGS analyses, one to four individuals per population were analysed. For 25 populations, ITS was sequenced and one sequence was included from GenBank, for 23 RGS was measured and 48 specimens from 39 populations were included in morphometric analyses (for details see Online Resource 1 and Fig. 1 as well as the text given as follows). In four cases, three specimens and in one case two specimens from the same collection were studied morphologically, which was important especially in the case of E. orjeni, for which only two collections were available.

DNA extraction, ITS sequencing and analyses of sequence data
Extractions of total genomic DNA and ITS sequencing were performed as described by Frajman and Schönswetter (2011), with the exception that sequencing was carried out at Eurofins Genomics (Ebersberg, Germany). In all cases, only a single individual per population was sequenced, exception being population 42 from Gramos, from which four individuals were sequenced, one morphologically corresponding to E. heldreichii and three morphological intermediates between E. amygdaloides and E. heldreichii. Contigs were assembled, edited and sequences aligned using Geneious Pro 5.5.9 (Kearse et al. 2012). Base polymorphisms were coded using NC-IUPAC ambiguity codes. GenBank numbers are given in Online Resource 1, and the alignment is available in Online Resource 2. Outgroup species were selected from all major clades resolved in E. sect. Patellares (Riina et al. 2013), including E. characias s.l., which occurs in the area of our focal taxa, E. heldreichii and E. orjeni, as well as the more distant E. erubescens and E. macroceras; the latter was used for rooting.
Maximum parsimony (MP) and MP bootstrap (MPB) analyses were performed using PAUP 4.0b10 (Swofford 2002). The most parsimonious trees were searched for heuristically with 100 replicates of random sequence addition, TBR swapping, and MulTrees on. All characters were equally weighted and unordered. The data set was bootstrapped using full heuristics, 1000 replicates, TBR branch swapping, MulTrees option off, and random addition sequence with five replicates. Bayesian analyses were performed using MrBayes 3.2.1 (Ronquist et al. 2012) applying the GTR substitution model proposed by the Akaike information criterion implemented in MrAIC.pl 1.4 (Nylander 2004). Values for all parameters, such as the shape of the gamma distribution, were estimated during the analyses. The settings for the Metropolis-coupled Markov chain Monte Carlo process included four runs with four chains each (three heated ones using the default heating scheme), run simultaneously for 10,000,000 generations each, sampling trees every 1000th generation using default priors. The posterior probabilities (PP) of the phylogeny and its branches were determined from the combined set of trees, discarding the first 1001 trees of each run as burn-in.

Relative genome size measurements
RGS was measured using flow cytometry (FCM) as described by Suda and Trávníček (2006). Nuclei of silica gel-dried material of the reference standard (Bellis perennis L., 2C = 3.38 pg; Schönswetter et al. 2007) and of the target Euphorbia sample were stained using 4′,6-diamidino-2-phenylindole (DAPI). The RGS was estimated for one to four individuals per population for 23 populations, 18 of E. amygdaloides, three of E. heldreichii and two of E. orjeni (see Online Resource 1). A CyFlow space flow cytometer (Partec, GmbH, Münster, Germany) was used to record the relative fluorescence of 3000 nuclei, and FloMax software (Partec) was used to evaluate histograms and to calculate coefficients of variation (CVs) of the standard and sample peaks. The RGS was calculated as the ratio between the values of the mean relative fluorescence of the sample and the standard. Mean value and standard deviation were calculated for each population.

Morphometric analyses
Twenty-eight individuals from 25 populations of E. amygdaloides, 14 individuals from 12 populations of E. heldreichii and six individuals from two populations of E. orjeni were analysed morphometrically (see Online Resource 1 for details).
Thirty-five metric characters were measured or scored; based on them, 12 ratio characters were calculated (Table 1). Leaf characters were measured on mid-stem leaves and raylet leaves from the basal bifurcation of rays. Plant height, stem length and width, stem leaf characters, number and length of axillary and terminal rays as well as number of branchings of terminal rays were measured manually. All the other characters (ray and raylet leaf, cyathium, fruit, seed and trichome characters) were measured on microscopic images taken with a stereomicroscope Olympus SZX9 using the Olympus image analysis software analySIS pro. Since fruits and seeds were developed only in a limited number of specimens (see Online Resource 1), we performed statistical analyses separately for this set of characters. Since fruits and seeds were not available in two and four specimens (two of E. heldreichii and four of E. amygdaloides), respectively, we replaced their character values in the data matrix with mean values of the respective taxon.
Statistical analyses were performed using Statistica vers. 5.1. (StatSoft 1996). Correlation among metric characters was tested employing Pearson and Spearman correlation coefficients but as correlation coefficients did not exceed 0.95 for any character pair, all measured characters were retained for further analyses. Box plot diagrams were produced for all characters, in order to visualise and show the variation among the three species. After standardisation to zero mean and one unit variance, principal component analysis (PCA) was performed. As a Tukey's HSD post-hoc test showed no discriminatory importance of fruit and seed characters, as well as of 17 other characters (LYSL, BLASL, LW, TrN, RyN, RyLL, RyLMW, RLL, CL, IL, IW, CGAL, including five ratios LMW/LL, RLL/RLW, RLMW/RLL, IL/IW, GL/GW), we excluded them from canonical discriminant analysis (DA). This analysis was applied to explore morphological differentiation between E. amygdaloides, E. heldreichii and E. orjeni and to clarify the relative importance of characters as discriminators between them. Values presented in the species descriptions and in the identification key correspond to the 25 and 75 percentiles, supplemented by extreme values in parentheses.

ITS phylogeny
The ITS alignment was 693 characters long, which was also the length of all ingroup sequences. Sixteen characters (2.3%) were parsimony informative. Homoplasy index was 0.05 (0.06 after exclusion of uninformative characters), and the retention index was 0.99. Totally, 5301 most parsimonious trees were found and their score was 20. Bayesian and maximum parsimony reconstructions resulted in congruent topologies (Fig. 2). One main clade excluding E. macroceras, which was used for rooting, was resolved (posterior probabilities, pp 1, MPB 98%), including in a polytomy the outgroups E. erubescens and E. characias, as well as two clades of ingroup taxa. One (pp 1, MPB 82%) included all accessions of E. amygdaloides, whereas the other (pp 0.75, MPB 65%) included one clade with E. heldreichii (pp 1, MPB 95%) and one with E. orjeni (pp 0.99, MPB 66%).

Morphometry
The morphological character states are presented in Online Resource 3. Box plot diagrams of the most important differential characters are shown in Online Resource 4. The PCA scatter plot (first three axes explaining 21.5, 11.7 and 10.0% of the total variation; Fig. 4a) based on vegetative and cyathium characters showed a separation trend with slight overlap between E. orjeni and E. amygdaloides/E. heldreichii along the first axis and between E. heldreichii and E. amygdaloides/E. orjeni along the second axis. The characters contributing most to the first separation were habitus characters (Plant height, Stem width, Number of branchings of terminal rays, Number of branchings of the longest Ratio distance from the base to the widest part of the fruit/fruit length FMW/FL axillary ray, Length of the longest axillary ray) and those to the second separation were Number of axillary ray whorls, Width of a ray leaf, Length of a raylet leaf, Width of a raylet leaf, Length of cyathial involucre. Accordingly, the DA of the same characters (Fig. 4b) showed a clear separation of all three taxa. Euphorbia heldreichii was separated from the other taxa along the first axis, and the characters contributing most to this separation were Number of axillary ray whorls, Ratio of Distance from the base to the widest part of a ray leaf/Length of a ray leaf, Number of branchings of terminal rays and Width of a ray leaf. Euphorbia orjeni was separated from the other two taxa along the second axis, and the most important characters contributing to this separation were Stem length, Plant height and Number of branchings of the longest axillary ray. For fruit and seed characters, both PCA and DA showed a strong overlap among the three taxa (not shown).

Discussion
Many Tertiary species found a shelter in the Mediterranean Basin that is recognised as one of the most important refugia for species of this ancient flora (Thompson et al. 2005) and as a biodiversity hotspot with an exceptional role in the preservation of unique species and genetic diversity (Myers et al. 2000;Nieto Feliner 2014   in polytomy with E. characias and E. amygdaloides. The sister relationship between E. heildreichii and E. orjeni is supported by only one synapomorphic base in ITS sequences (alignment position 625) and their divergence, supported by several differences in ITS sequences, happened later, possibly in the Pliocene. All three species survived the Pleistocene glaciations in the Balkan Peninsula but reacted differently to Holocene warming.
Euphorbia amygdaloides is distributed throughout Europe, North Africa and Western Asia (Fig. 5), suggesting that it might have had multiple Pleistocene refugia, from which it extended its range considerably in parallel with the spread of forests in the Holocene, similar to other understory forest species in this part of Europe (Kramp et al. 2009;Willner et al. 2009;Slovák et al. 2012;Rešetnik et al. 2016).
On the other hand, E. heldreichii, the most thermophilous of all three species, remained confined to the southern Balkan Peninsula, similar to Aesculus hippocastanum and Forsythia europaea (Xiang et al. 1998;Ha et al. 2018). In its northern distribution area in Albania and northern Greece, its range overlaps with that of E. amygdaloides (Aldén 1986;Barina 2017;Fig. 5). Both species can occur in close vicinity, e.g. in the region of Mt. Gramos in Greece (Fig. 1), and Aldén (1986) reported that intermediate individuals occur in this area. However, forms of E. heldreichii with less pronounced whorls of axillary rays from Mt. Gramos (population 42; three individuals sequenced), thus resembling E. amygdaloides, had the same ITS sequence as co-occurring typical E. heldreichii (one individual sequenced). This suggests that they are probably not hybrids between E. amygdaloides and E. heldreichii, but just atypical forms of the latter and that the two species do not hybridise.
Euphorbia orjeni (Fig. 6), the most enigmatic of the three species, is the most restricted. It is only known from its type locality below Mt. Orjen in Montenegro and from one other Leaf material was sampled from this specimen deposited at BEOU for DNA extraction due to its monstrous appearance and several prolonged axillary rays, and its ITS sequence revealed that it belongs to E. orjeni. It is thus possible that E. orjeni has a larger distribution in this part of the Balkan Peninsula (Montenegro, Serbia, and possibly adjacent Bosnia and Herzegovina). Despite hosting a big collection of specimens from the Balkan Peninsula, none of the herbarium specimens of E. amygdaloides from the Balkans deposited at W and WU morphologically corresponded to E. orjeni, indicating that the species is rare. However, due to its resemblance to E. amygdaloides, which is a common species and thus not attractive for herbarium collectors, and its in general bigger and more branched habit, which is impractical for preparing herbarium specimens, it has likely remained overlooked and undersampled. Future, more focused search for E. orjeni will reveal where the species is actually distributed.
Our morphometric survey has shown that the three species are also morphologically distinct and confirmed the discriminatory power of characters that were traditionally deemed to distinguish between E. amygdaloides and E. heldreichii, i.e. the presence of whorled leaves and axillary rays on the stem. Similarly, our study confirmed some of the characters that were considered diagnostic for E. orjeni in relation to E. amygdaloides by Beck Mannagetta (1920). However, since Beck Mannagetta (1920) inspected only two flowering shoots lacking the basal parts, some characters stated in the protologue show more variation than suggested. For instance, the nectar glands are not always green but can also be yellowish or light brown, and their appendages/horns are not significantly longer than in the other two species and are only sometimes dilated and emarginate at the tip. In addition, the plants are not always glabrous, but can sometimes be sparsely pubescent. Its leaves are somewhat thicker and coriaceous compared to E. amygdaloides, and thus similar in consistence to E. heldreichii. Finally, the plants did not have a notably unpleasant odour as stated by Beck Mannagetta (1920: "odore fastidioso"), probably based on the Veselý's note on the herbarium label of the type specimen ("sehr übel riechend"), and were mostly only up to 0.5 m high, and not "c. 1 m", as stated by Beck Mannagetta (1920), also based on herbarium label ("1/2 mannshoch"-"half-man high"), as he had not seen the entire plant ("partes inferiores non vidi").
The RGS data of the three species were not significantly different, as the RGS values of E. heldreichii and E. orjeni were within the range of those of E. amygdaloides, which exhibited considerable variation in RGS among populations and a certain trend of increasing RGS from the east towards the west. The RGS of E. heldreichii was lower than that of E. orjeni and was in the range

Taxonomic treatment
Based on our phylogenetic and morphological data, we provide the following taxonomic treatment, recognising E. amygdaloides, E. heldreichii and E. orjeni as independent species.
Conservation status: Following the criterion D of the IUCN (2012) for vulnerable species, i.e. population size estimated to number fewer than 1000 mature individuals with a very restricted area of occupancy, we deem E. orjeni vulnerable (VU) based on the current knowledge of its distribution.