Morphometrics as a robust tool for disambiguation in plant taxonomy: the case of Lactuca livida, a commonly accepted but never delimited taxon

Abstract An important goal of taxonomy is to clarify the identity of ambiguous taxa. Since lineage divergence usually involves ecological shifts that are associated with plant morphology, we propose that searching for fixed, non-overlapping morphological characters and specific associations of features should be the first task in disambiguation. We applied this idea in an analysis of the taxonomic identity of Lactuca livida, an Iberian endemic relative of L. virosa that, despite an extremely imprecise delimitation, is usually recognized in standard floras and germplasm banks. We analysed 24 possible diagnostic characters across the Iberian Peninsula, drawn from related taxonomic literature. No discontinuities in frequency variability distribution characters were found, even in two bimodal quantitative characters: leaf lobation and number of florets per capitulum. There were no notable patterns of association among characters, and the PCA/PCoA score plots did not show any distinctive groupings. Leaf lobation followed a significant geographic pattern, but there was no effective segregation of leaf shapes. We conclude that the variability found is symplesiomorphic, as it is present throughout the Serriola group of Lactuca, to which L. virosa belongs. The analyses, together with previous biological knowledge, indicates that Lactuca livida should be considered a synonym of L. virosa.


Introduction
Accurate species recognition and delimitation is one of the main goals in the practice of taxonomy. Over the last century, this practice has become progressively more complex in the framework of an intense debate about the nature of species, which has resulted in about 30 different species concepts (de Queiroz 2007;Wilkins 2009;Zachos 2016), each giving rise to a particular taxonomic methodology (Luckow 1995;Marshall 2003, 2004;Agapow and Sluys 2005; but see McDade 1995). Aside from the old typological approach, the only species concept that clearly promotes morphometric analysis is the phenetic species concept, which is tested by analyzing the mathematical grouping pattern of the values of diverse morphological characters among a set of individuals or taxonomic units (Michener and Sokal 1957;Michener 1970;Sneath and Sokal 1973;Sokal 1986). The recognition of morphological entities as the result of an evolutionary biological process connects morphometrics to phylogenetic theory and the phylogenetic species concept (Nelson and Platnick 1981;Cracraft 1983;de Queiroz and Donoghue 1988;Nixon and Wheeler 1990;Crisp and Weston 1993;Baum and Donoghue 1995;Luckow 1995;Henderson 2006), thus, a deeper exploration of possible collaborative approaches between morphometrics and other methodologies would be desirable (Wiens 2007).
Currently, the general tendency in taxonomy is to support species delimitation in the case of independently evolving lineages (Fujita et al. 2012;Cutter 2013;Conix 2018; but see Freudenstein et al. 2017). At the same time, it has become an increasingly common practice to seek different lines of evidence of taxa differentiation that correspond with phylogenetic outcomes in order to promote taxonomic stability; this has been called integrative taxonomy (Dayrat 2005;Schlick-Steiner et al. 2010;Conix 2018;Dantas-Torres 2018). In this context, morphometric analysis and phenetic characterization are not considered the most robust evidence for taxonomic decisions because they are not involved on the ontological debate about the nature of species (Valcárcel and Vargas 2010). In addition, they cannot detect processes of homoplasy, either by parallel evolution or by evolutionary convergence (Arendt and Reznick 2008;Pearce 2012;Mejías et al. 2018).
In practice, however, morphological analysis remains the most frequently used methodology in taxonomy (Stevens 1991(Stevens , 2000McDade 1995;Knapp et al. 2005;Wiens 2007;Stuessy 2009; Dantas-Torres 2018) and has a considerably high level of success compared to other methods (Schlick-Steiner et al. 2010). Morphological analysis can be used alone (e.g., Henderson 2004;Valcárcel and Vargas 2010;Jiménez-Mejías et al. 2014;Morales et al. 2014) or in combination with other disciplines (e.g., Koffi et al. 2010;Ronikier and Zalewska-Gałosz 2014;Prata et al. 2018;Sokoloff et al. 2019). In general, morphology is essential in the practice of taxonomy for several reasons: i) morphologic (phenetic) discontinuities provides well-founded evidence of reproductive isolation (Stace 1989;Sites and Marshall 2004) and population genotypic aggregation (Mallet 1995(Mallet , 2020, since it is commonly assumed that morphological differentiation has a genetic basis (Wiens 2007;Valcárcel and Vargas 2010); ii) it is the simplest, cheapest and most direct method in taxonomy (Stace 1989); and iii) in combination with ecological approaches, it is quite helpful in understanding patterns of phenotypic variation and processes that lead to speciation (e.g., Turchetto et al. 2014;Freudenstein et al. 2017;Robbiati et al. 2017;Vogel Ely et al. 2018). In addition, morphology is the most useful way to identify diagnostic characters for species (Stace 1989;Valcárcel and Vargas 2010;Sokoloff et al. 2019), which facilitates the role of taxonomy as the supplier of basic units of study in ecology, biogeography, evolutionary biology and conservation (Brown et al. 1996;Blackburn and Gaston 1998;Barraclough and Nee 2001;Mace 2004;Balakrishnan 2005). Thus, from a modern taxonomic perspective, a deep and updated study of morphology is essential for both the practical delimitation of taxa and detection of the phylogenetic patterns that give rise to taxonomic diversification (Wiens 2007).
A particularly challenging goal in common taxonomic practice is to clarify the identity of numerous persistent problem taxa, some of which were described long ago and have never been precisely delimited (Patterson et al. 2006;Pometti et al. 2007;Potapova and Hamilton 2007;Kougioumoutzis et al. 2016;Vogel Ely et al. 2018;Dauncey et al. 2019;Duan et al. 2019). From our perspective, the first basic analysis to be done is a detailed morphological prospection to detect characters or particular associations of traits that may have gone unnoticed and that provide evidence of differentiation processes (Mallet 1995;Wiens 2007;Valcárcel and Vargas 2010). Lactuca livida Boiss. & Reut., a close relative of the common weed L. virosa L., is one of these troublesome taxa in the Iberian flora (Lindqvist 1960b;Feráková 1977;Mejías 1993Mejías , 2017Lebeda et al. 2001. It was separated from L. virosa by E. Boissier and G.F. Reuter (in Boissier 1845), mainly because of particularities in the coloration, leaf morphology and number of florets per flower head (Table 1, Figure 1) of some materials collected in Montes de Toledo (central Iberian Peninsula). However, morphological boundaries between L. virosa and L. livida have remained unclear and imprecise ever since, giving rise to some taxonomic inconsistencies that are of particular interest for the compilation of available genetic resources for lettuce breeding Sretenović Rajičić et al. 2008;Doležalová et al. 2003Doležalová et al. , 2004. Feráková (1976Feráková ( , 1977, in the revision of the genus Lactuca L. for Flora Europaea, recognized the validity of the two taxa, but pointed out the strong similarities between them and kept the specific taxonomic rank for both. Some other authors have also used this taxonomic treatment in floras of narrower geographic scopes (e.g., García Rollán 2005). On the contrary, Velasco (1981) proposed that L. livida should be subordinated to L. virosa L. at the subspecific level, so he coined the combination L. virosa subsp. livida (Boiss. & Reut.) Ladero & Velasco. Moreover, he suggested that this taxon should be present in some other mountain areas of the Iberian Peninsula. Cueto (1985, 2009) used this new combination L. virosa subsp. livida and indicated some new localities in Southeastern Spanish ranges, which reinforced the support for Velasco's proposals. Subsequently, it was reported that the two species have the same number of chromosomes and show no notable differences in chromosome morphology (Mejías 1993), although these observations conflict with some variations in DNA content and puzzling genetic diversity patterns detected in samples from germplasm collections (Doležalová et al. 2002, Sretenović Rajičić et al. 2008. No differences in the reproductive system between the two proposed species have been detected (Mejías 1994).
Despite the reported botanical and agronomic interest in disentangling the taxonomic ambiguity between Lactuca virosa and L. livida, as far as we know, no phenetic or phylogenetic study has been performed. The first step toward the clarification of the validity of Lactuca livida is a detailed morphological analysis elucidating whether any phenetic group showing the characters attributed to this taxon can be separated within L. virosa s.l. in the geographic range of the Iberian Peninsula. Since no karyological and/or reproductive differentiation has been detected within L. virosa s.l. in the area (Mejías 1993(Mejías , 1994, it is reasonable to assume that a hypothetic speciation process might be mediated by particular ecological constraints (Schluter 2001;Nosil 2012;Favre et al. 2017). In line with this, our present goals are: (1) to record the morphological variability of L. virosa s.l. in the Iberian Peninsula in order to detect non-overlapping morphological characters and/or particular association of features, with particular emphasis on characters considered to be  Velasco (1981), ( 5 ) cueto (1985, 2009 6-10 mm, blackish, 5-ribbed 3 up to 7 mm 3 , blackish 1, 3 , 5-ribbed 1 or 5-6-ribbed 3 Achene beak as long as the body 3 as long as the body or shorter 3 Figure 1. type material of Lactuca livida Boiss. & reut. and details of morphological variability of L. virosa l. sensu lato on the iberian Peninsula with their scores. a. lectotype of Lactuca livida (g00300061, specimen on the left side) and one isotype (on the right side), assigned by Burdet et al. (1983). b. Voucher specimen (Vianos-tortas, seV285532) showing basal leaves with long petiole, bottom stem colour = 4 and lobation coefficient = 0.15. c. Voucher specimen (torres, seV286663) with bottom stem colour =1 and lobated, runcinated leaves (lobation coefficient = 0.83). d. midrib thorn abundance = 3, leaf colour = 1 (albergaria-a-Velha, seV286648). e. upper stem colour = 4, upper stem thorn abundance = 2 (sierra del caurel, seV565679). f. upper stem colour = 1, upper stem thorn abundance = 4 (santander, seV ma680349). g. midrib thorn abundance = 4, leaf nerve thorn abundance = 4, leaf colour = 3 (serra de são mamede, seV286641). h. ligule colour = 2, inflorescence calyx colour = 2 (carcabuey, seV286659). i. ligule colour = 4, inflorescence calyx colour = 4 (rozas de Puerto real, ma 556961). j. mature infructescence with achenes to be dispersed (carcabuey, seV286659). k, l. Variability in achene body shape (k: san nicolás del Puerto, seV286650; l: san Pablo de los montes, seV285540). m. achene ribs n = 8 (san Pablo de los montes, seV285540). discriminative of L. livida in the literature (Marhold 2011); 2) to analyse geographical and ecological patterns associated with morphological variability to bring to light possible selective forces promoting differentiation processes; 3) to discuss taxonomic implications of the morphometric analysis performed and other previous biological knowledge of the group, according to relevant alternative species concepts and current taxonomic tendencies; and 4) to exemplify the role of morphometrics in taxa disambiguation and its impact on plant taxonomy applications.

Plant materials
Our study was mainly based on measurements taken from 164 herbarium specimens (Online Supplement Appendix S1) collected throughout the Iberian Peninsula (155) and southern France (9). Most commonly, we used one or two individuals per population, but in some cases (e.g., population of San Pablo de los Montes, Toledo) we studied up to five, often collected in several years. Specimens are currently kept in the herbaria ARAN, BC, GR, HuAL, JA, LEB, LISE, LISu, MA, SALA, SALAF, SEV and VAL. In addition, we analysed the variation in the number of florets per capitulum in eight populations (see Table 7). For that, we randomly collected developed flower head buds in the field or in the experimental garden, which were preserved in a 70% ethanol solution to be subsequently counted in the laboratory. We performed two types of sampling in each population: variation in a single individual, for which we collected 15 to 51 flower heads from one randomly selected individual, and variation throughout the population by random selection of 10 to 20 individuals per population, from which we collected one head each.

Morphological characters
According to the literature reviewed, twenty-four characters were considered to be important to discriminate L. virosa L. from L. livida Bois & Reut.; twelve of them were quantitative and the remaining twelve qualitative (Table 2). In each specimen, we always selected well-developed and well-preserved plant organs for the study. Measurements were taken with a Mitutoyo digital calliper (accuracy ± 0.01 mm). A Leica MC170 HD stereoscopic microscope was used to determine the values of micro-characters.

Geographical information
The geographic coordinates and altitudes of collection for the specimens studied were drawn from the label attached to the exsiccatum (herbarium specimen). When coordinate data were missing, they were estimated by searching for the location name given on the label using GoogleEarth 7.1.3.22.3 (https://google_earth.es.downloadastro.com/old_versions/). All geographic coordinates were transformed to decimal degrees, rounded to the nearest tenth of a degree.

Statistical analyses
Statistical analyses were run in R version 3.5.1. The following packages were employed: base package (R Core Team 2018), readxl ( (Huang et al. 2014). Table 2. taxonomic characters used in the study and the criteria to measure them.

Quantitative Characters Criteria (Units)
Number of florets number of florets in a randomly selected capitulum (n) Achene length achene body length of a randomly selected achene (mm) Achene width achene's body biggest width (mm) Achene ribs number of ribs present in the achene's body (n)

Peak length
Peak's length of the achene (mm)

Pappus length
Pappus' hairs longest length (mm) Basal petiole length length of the petiole in a basal leaf (mm) Stem thorn length length of the longest thorn in the stem (mm) Leaf thorn length length of the longest thorn in the leaf (mm) Leaf largest width maximum width near the centre of the biggest non-basal leaf (mm) Leaf smallest width minimum width near the centre of the biggest non-basal leaf (mm) Lobation coefficient 1 -(leaf smallest width ÷ leaf highest width)

Qualitative Characters Criteria (Range)
Achene colour colour of a mature achene ( Normality of distribution of characters/traits was tested by the Anderson-Darling test when the sample size was greater than thirty and by the Shapiro-Wilk test in the case of lower sample sizes. For variables suspected to be bimodal, the Hartigans' Dip Test for unimodality was computed and the modes, antimodes, amplitudes (distance between modes), bimodality ratios and P-values of the contrasts were calculated.
In order to test the association among quantitative and qualitative characters, quantitative data were grouped by the scores of the specimens on the qualitative scales (groups 1 ~ 4). Given that we only recorded five states for the discontinuous quantitative character achene ribs, it was also used as a grouping criterion for the remaining quantitative characters, similar to the qualitative characters. The means of the resulting groups were tested for equality to detect potential distinctive groups by running one-factor ANOVA and Kruskal-Wallis analyses. Tukey-HSD and Pairwise Wilcoxon tests were applied as post-hoc analyses when necessary, under Benjamini-Hochberg's procedure for correction of false discovery rate (FDR) for type I errors. We explored if there were distinctive groups by performing clustering and ordination multivariate analyses. We excluded the variables basal petiole length, bottom stem thorn abundance and bottom stem colour from all analyses because they had many missing values and dramatically decreased sample size. We also removed lobation coefficient from the analyses, since it is calculated from two existing variables. No characters reached a high enough correlation value to justify their exclusion from the analyses (|r| or |ρ| < 0.90; Table 5). Therefore, we analysed the euclidean distance matrix of ten quantitative characters using hierarchical cluster analysis (uPGMA), and PCA analysis. Hierarchical cluster was tested using a significance analysis of clustering. PCA was evaluated through Kaiser-Meyer-Olkin statistics (KMO) and Bartlett's Test of Sphericity; consistently, the variable achene ribs was excluded to preserve the adequacy of PCA sampling according to the KMO values (see Results). We analysed the Gower distance matrix of ten quantitative and ten qualitative characters using a PCoA analysis without correcting negative character loadings. For graphical representation, we scaled the eigenvectors to the square root of the corresponding eigenvalues. This way, we could preserve the original distance among plant specimens (Gower 1966). We also run a hierarchical cluster analysis based on the distance matrix used in the PCoA. Correlations between altitude and quantitative characters were checked by Spearman's correlation tests. For the two traits that had bimodal distributions (number of florets per capitulum, and leaf lobation coefficient), we determined whether the trait values followed geographic patterns using Mantel tests with one thousand permutations. Similarly, we combined the remaining characters to generate five variables to test other possible geographic patterns (Online Supplement  Table S4).
Comparisons between intraindividual and intrapopulation variabilities for the number of florets per flower head were addressed using a Wilcoxon rank test for paired range amplitude measurements in each population in order to detect possible variation of the character within the populations.

Frequency distribution of characters
Normality tests showed that seven out of the twelve quantitative diagnostic characters were not normally distributed (Table 3). After plotting them as histograms, (Online Supplement Figure S1), the following characters were suspected to be bimodal: number of florets per capitulum, achene length and leaf lobation coefficient. Hartigans' Dip Test scores for unimodality (Online Supplement Table S1) indicated that number of florets (P < 0.001) and leaf lobation coefficient (P = 0.013) were bimodal variables with frequency distributions originated by overlapping two normal distributions ( Figure 2). For number of florets, n = 15 (15.19) and n = 19 (18.70) were the modes and n = 17 (17.68) was the antimode, while for leaf lobation coefficient the modes were 0.22 and 0.84 and the antimode was 0.56. A short amplitude (A = 0.054) for florets per flower head showed high overlap of the normal distributions, while for leaf lobation coefficient (A = 0.473) both distributions appeared to be well-defined. Achene length was found to be unimodal (P = 0.963).
ANOVA and Kruskal-Wallis tests indicated several differences in quantitative characters among groups of specimens arranged according to qualitative character states (Table 4), but post-hoc tests failed to find any distinctive groups. The analyses showed that the means differed, especially for thorn-related variables. Discrepancies in thorn abundances and lengths were expected since there was a group of spineless plants, which consequently exhibited a mean spine length equal to zero (Online Supplement Figure S1-g). Nevertheless, post-hoc results also showed that plants with a high density of thorns tended to possess longer thorns on both the stems and leaves. In the cases of the other significant ANOVA, we interpret them as false positives, since the post-hoc tests often did not show differences between the extreme states, but rather patterns in which intermediate states were found to be significantly different from the extremes, which do not have a clear interpretation (Online Supplement Table S2). We also found no association between the number of achene ribs (a quantitative discontinuous character with only five states) and the remaining quantitative variables. Table 3. exploratory analysis of quantitative characters. range, mean (μ) and standard deviation (σ) are displayed along with P-values of the normality tests that were performed on each one: shapiro-wilk (Basal petiole length) and anderson-Darling. Variables with a * were shown to be bimodal by later analyses.

Cluster and ordination analysis
Correlation coefficient contrasts revealed several significant correlations among quantitative characters (Table 5). Therefore, ten variables were chosen for inclusion in the hierarchical cluster analyses (uPGMA, Table 6). We removed missing values (N = 66), and calculated a euclidean distance matrix based on the remaining data. The hierarchical cluster dendrogram showed one branch of five plant specimens, coming from very disparate geographic locations in the study area, separated from the others, which in turn were subdivided in two even branches (Online Supplement Figure S3). Nevertheless, significance analysis of cluster revealed there is no grouping pattern (P = 0.567). Initially, the KMO test score (KMO = 0.488) showed an unacceptable adequacy to run the PCA; however, after removing the variable achene ribs (KMO = 0.534) it was sufficient (Kaiser 1974). Bartlett's Test of Sphericity (P < 0.001) indicated that the data (nine variables) were suitable for PCA. Variables related to achenes had high weights in component 1, while thorn lengths exhibited high weights in component 2; surprisingly, the weight of flowers per capitulum was moderate in component 1 and low in component 2 ( Table 6). The PCA ordination diagram ( Figure 3A) showed that the plant specimens formed a continuous spectrum, with more individuals at the extremes. We selected twenty-two quantitative and qualitative characters for the PCoA analysis. We excluded missing values (N = 53), and calculated a Gower distance matrix based on the remaining data. The PCoA plot ( Figure 3B) showed that the plant specimens formed a continuous spectrum similar to that observed in the PCA analysis. The hierarchical cluster based on the Gower distances did not reveal any grouping pattern (Online Supplement Figure S4).

Geographical patterns
No significant correlations were found between altitude and the eleven quantitative traits recorded across the Iberian Peninsula (Online Supplement Table S3). The Mantel test showed that leaf lobation coefficient values followed a geographic pattern (P = 0.008), while the geographic distributions of the number of florets per capitulum (P = 0.575) and the remaining variables tested (P = 0.187 − 0.768) were random (Online Supplement Table S4). Mapping the leaf lobation coefficient (Figure 4) revealed that plants tended to have a more pronounced leaf lobation (coefficient > 0.5) in the western Iberian Peninsula, and plants located on the east side of the peninsula and in southern France exhibited mild lobation or no lobation (coefficient ≤ 0.5).

Variability in the number of florets per flower head within populations
Variation in the number of florets per capitulum within individuals and populations were analysed in eight different locations (Table 7). Variability at both individual and population levels was found to be quite high, given that the general range of variation found throughout the Iberian Peninsula was (9)11-23 (Table 3, Figure 2). No significant differences in range amplitude (P = 0.235) or variation coefficient (P = 1.000) were detected for this variable at either level of analysis. Surprisingly, the antimode value (n = 17) from the bimodal distribution obtained in the general sampling has been found to be within the ranges calculated for the intraindividual and the intrapopulation variability, with the exception of the samples from Sierra Nevada, where the largest number of florets was 15.

Discussion
In general, we found wide and continuous variability in the characters commonly used in the literature as criteria to distinguish between Lactuca livida Boiss. & Reut. and L. virosa L. throughout the Iberian Peninsula (Tables 1 and 2). It is also notable that the traits considered to be specific to L. livida were common throughout the area (see Online Supplement Figures S1 and S2). However, despite including material from the type locality of L. livida and from other areas where it was expected to be present (Velasco 1981), we did not detect non-overlapping or diagnostic morphological traits, nor did the PCA/PCoA show distinctive grouping of specimens (Figure 3). From a phenetic standpoint, our results clearly show the difficulty in detecting morphological grouping within the plants studied, which is not surprising in a group as persistently taxonomically problematic (Wiens 2007) as L. livida. From a cladistic perspective, our results reveal the apparent absence of derived or autapomorphic diagnosable traits in the area, which indicates that any possible speciation processes must be due to a particular combination of plesiomorphic features (Cracraft 1983;Donoghue 1985;de Queiroz and Donoghue 1988;Olmstead 1995). Considering this rationale, and that clear diagnosability is not actually a requirement of most species concepts (Wiley 1978;Nixon and Wheeler 1990; but see Cracraft 1983), we discuss the variability patterns recorded and possible associations among particular traits (Nixon and Wheeler 1990;Davis and Nixon 1992;Mallet 1995;Grant and Grant 2006;Henderson 2006;Valcárcel and Vargas 2010) in order to detect possible differentiations of taxonomic relevance within L. virosa s.l. on the Iberian Peninsula. Moreover, we integrate morphometrics with geographic and altitudinal analyses to assess possible ecological patterns of variability denoting morphological ecotypic differentiation (Clausen 1967;Brown et al. 1996;Gianoli et al. 2004 Table 4. P-values of anoVa and Kruskal-wallis test results for each quantitative character, grouped by each qualitative character. several differences in means were found, but post-hoc analyses failed to find any distinctive groups except for stem-related characters. however, the differences in these characters were caused by one of the groups being formed by thornless plants or by the positive association between longer thorns and higher thorn density. non-normal quantitative characters are marked with *. It is surprising that among the twelve quantitative characters analyzed, all but one (stem thorn length) showed continuous frequency distributions (Figure 2, Online Supplement Figure S1). Moreover, several of them conformed to a normal distribution (Table 3) and only two showed bimodal distributions (Figure 2, Online Supplement Table S1). That denotes a limited value of this set of characters as indicators of a hypothetical differentiation process and for taxonomic diagnosis. The discontinuous distribution of stem thorn length reflects a notable abundance of specimens with no thorns on the stem, and the general association between the spination traits shows a clear tendency of plants with higher thorn density to have longer thorns both on the stem and on the leaves (Tables 4 and 5). We did not detect any general association of spination traits to other quantitative plant features with the exception of the normal variability of achene width (Table 5). In addition, no clear geographical or ecological (altitude) pattern was identified (Online  Table 7. intra-individual and intra-population variabilities of the quantitative character Number of florets illustrated by the range, the mean (µ) and the standard deviation (σ). a wilcoxon test for range breadth showed that the degree of intra-individual variation did not differ from the intra-population variation (see P-values). samples were collected in the wild, except those with an asterisk (*), which were obtained from plants growing in the experimental garden. e: spain, P: Portugal.  Supplement Tables S3 and S4). Therefore, despite their prevalence in the literature, it seems unlikely that thorns are of taxonomic relevance here, and no autapomorphic traits can be identified among the characters involved in spination, at least on the Iberian Peninsula. Lindqvist (1960b) had already pointed out high variability of spination in L. virosa and reported a similar situation in the L. serriola-L. sativa group, which led to his proposal of a genetic system of one major gene with four alleles and some modifiers as the mechanism of inheritance (Lindqvist 1960c). We hypothesize that this variability in spination is a symplesiomorphic feature and that a genetic system similar to that of L. serriola and L. sativa could be acting in L. virosa. None of the five quantitative fruit size features analyzed provide evidence of morphological differentiation. All showed unimodal distributions, three of them conforming to normality (Table 3, Online Supplement Table S1), and they had limited correlations to other groups of traits (Tables 4 and  5). An exception was that the clearly normally distributed variables achene width and pappus length are associated to stem thorn length (as discussed above) and petiole length of basal leaves, respectively. However, these associations do not allow the identification of a diagnosable or recognizable group of plants. In addition, the absence of clear altitudinal and geographical patterns of variability (Online Supplement Tables S3 and S4) in these and the remaining unimodal quantitative characters does not provide evidence of any ecologically mediated process of morphological diversification.

Infloresce-nce type
Of the two variables that conform to bimodal distributions, leaf lobation coefficient proved to be of considerable interest. The frequency distribution found (Figure 2) indicates that two leaf shapes-deeply lobed and almost entire-are predominant in Lactuca virosa L. s.l. on the Iberian Peninsula. There is a significant geographic pattern of variability (Online Supplement Table S4), but the segregation between the types is quite diffuse; both shapes can be found throughout the range (Figure 4), and it is not uncommon to find both types in a single population (J.A. Mejías, pers. obs.). In several species, leaf shape has been shown to affect water supply trade-offs and thermoregulation, with lobed leaves considered advantageous in conditions of high irradiance (Vogel 1970;Sisó et al. 2001;Nicotra et al. 2011;Campitelli et al. 2013). Thus, we can hypothesize that the pattern detected reflects some kind of ecological differentiation. Since L. virosa s.l. colonizes lower-altitude areas exclusively towards the northern and northwestern parts of the Peninsula, and is present only in highland areas in the south, an association between altitude and leaf lobation was plausible. However, our statistical results do not confirm this hypothesis (Online Supplement Table S3) and the ecological and biogeographic significance of the pattern remains unclear. In addition, lobation variability is not correlated with any other character analyzed (Tables 4 and 5). A similar differentiation in leaf shapes is present in L. serriola L. (Mejías 2017), in which both leaf types can often be found in the same populations (Prince and Carter 1977;Lebeda et al. 2001;J.A. Mejías, pers. obs.). Based on an extensive breeding program, Lindqvist (1958) proposed that the inheritance of the two alternative shapes in L. serriola is determined by a single gene with two alleles, of which the lobed form is the dominant trait. It seems reasonable to hypothesize a similar genetic system for L. virosa and to propose that this variability is symplesiomorphic.
The number of florets per flower head also displays a bimodal distribution; the lower mode was 15, which matches indications for Lactuca virosa s. str. in the literature (Feráková 1976(Feráková , 1977, and the higher mode was 19, a less divergent value than the "c. 25" proposed by Boissier and Reuter (in Boissier 1845) and Feráková (l.c.) for L. livida (Table 1). Surprisingly, the character is not correlated to thorn length or leaf lobation coefficient (Table 4), two variables that are commonly used in the distinction of L. livida. Variation in the number of florets per head is not a well-studied character, so little is known about its genetic regulation. A few studies concerning number of ray florets in flower heads propose the involvement of major genes (Moritz and Kadereit 2001;Song et al. 2018). Our results seem to show a disomic inheritance mechanism based on one gene with dominance and several modifiers, which allows the concurrence in populations and individuals of highly variable numbers of florets ( Figure 2, Table 7). On the other hand, it is clear that a greater number of florets per flower head may increase the showiness of the blooms, but the range of variability observed seems unlikely to impact the pollination system in ways that would have evolutionary consequences. There were no trait-related geographical or altitudinal patterns that could shed light on any possible ecological significance (Online Supplement Tables S3 and S4), and again, differences in number of florets per flower head provided no clear evidence of a taxonomic differentiation process.
The variability of states in qualitative characters studied do not show any remarkable patterns of association to quantitative features, with the exception of thorn abundance, which, as expected, is clearly associated with thorn length (Table 4). According to the general comments in the literature (Table  1), we would have expected a clearly defined pattern of association among thorn traits, the number of florets per flower head and the incidence of purple coloration. However, we only found an association of purple stem coloration and leaf thorn length. Interestingly, Lactuca serriola and L. sativa also have variable patterns of anthocyanin pigmentation in a complex inheritance mechanism (Lindqvist 1960b(Lindqvist , 1960c) that depends on growth conditions (Brücková et al. 2016), which suggests that purple stem colour is another symplesiomorphic trait. Other scattered associations (Table 4) were somewhat unexpected. We wonder if they reflect ecological adaptations or are just random, since the grouping patterns were not consistent with any ecological and/or taxonomical differentiation between the extreme groups, but rather, differences between intermediates and extremes (see Online Supplement Table S2). We were not able to analyse life-form variability (Table 1) but, when present, voucher specimens commonly show a thickened taproot that indicates a clear trend towards biennial cycle (Feráková 1977), in contrast to the clearly annual relatives L. serriola L. and L. saligna L. According to our observations, wild individuals from Iberian populations seem to most commonly reach the flowering stage during the first year, so few vegetative rosettes remain at the end of the season. However, plants from these same populations grown in pots in our experimental garden usually behave as biennials. Time of flowering is an important adaptive trait that involves a range of environmental variables and, in general, there is no sharp distinction between annuals and biennials (Bouché et al. 2017;Amasino 2018). In this case, it is clear that plants and populations can modulate the life cycle expression. Therefore, although the geographic pattern of life cycle variability in this species probably merits further study, it seems unlikely to be of taxonomic interest.
The above discussion clearly shows that most traits assigned to Lactuca livida Boiss. & Reut. are within the range of current variability of L. virosa. The only discontinuity detected involves stem thorns, but the weak association with other traits does not allow the identification of any distinct group of specimens. Furthermore, the variability described is present in other species of the Serriola group (Lindqvist 1958, 1960a, 1960b, 1960c, Mejías 2017 to which L. virosa belongs to (Koopman et al. 1998(Koopman et al. , 2001Güzel et al. 2021). Hence, according to our analysis, no apomorphies can be identified, and putative specific attributes of L. livida can be mainly regarded as symplesiomorphies in the Serriola group. A notable exception is the number of florets per capitulum, which, as far as we know has not been analyzed before, but as we discussed earlier, does not allow us to differentiate diagnosable morphological groups or to propose the divergence of a new lineage. All of these conclusions are well reflected in the plots of the first two principal components of the PCA/PCoA analyses, where all specimens clump into a single group (Figure 3). Moreover, the specimen that we collected from the type locality (San Pablo de los Montes, Toledo) plotted close to the core of the group, although PCoA suggests that probably it does not show the most common combination of morphological features. Similarly, it shows central positions in hierarchical dendrograms (Online Supplement Figures S3 and S4). This indicates that this material does not represent a differentiated morphotype, or even a marginal combination of traits, as often occurs in widespread plants (e.g., Pometti et al. 2007;Vogel Ely et al. 2018). Indeed, it is not uncommon to detect plants showing one or more of the attributes considered to be specific to L. livida along with other traits of L. virosa s. str. (e.g., individuals bearing thorns and lobed, runcinated leaves but small number of florets per flower head). We were unable to analyze in precise detail the type material of L. livida designated by Burdet et al. (1983), but it also shows a quite common appearance of L. virosa of the Iberian Peninsula, with deeply lobed stem leaves, long entire basal leaves and spinulose stem and leaves (Figure 1). All of these findings align with an absence of support for the validation of L. livida under a phenetic species concept. The conclusion is the same under an original cladistic perspective (Cracraft 1983;Nixon and Wheeler 1990;Baum and Donoghue 1995;Olmstead 1995), since our results show that splitting into two taxonomic groups also seems unfeasible. The alternative possibility of identifying all peninsular plants as L. livida, including individuals from southern France, is totally unsuitable. It would mean that materials from the Iberian Peninsula and southern France constitute a different taxon from the rest of European plants, and no evidence has ever been reported to suggest this (Feráková 1976(Feráková , 1977. All of these considerations, however, seem to show the existence of some morphological gradients across Europe that could be analyzed from an ecological and geographical point of view. Plants from the type locality of L. livida and typical material of L. virosa both show a common chromosome number (n = 9, 2n = 18), with very similar chromosome morphology and karyotype characteristics (Lindqvist 1960a;Mejías 1993), unique in the Serriola group. This allows us to hypothesize that typical specimens of L. virosa and L. livida are able to easily cross and generate fertile offspring (Sokal and Crovello 1970;Levin 2002;Baltisberger and Hörandl 2016;Deanna et al. 2018). The strong karyological similarity is not consistent with the differences between the two species' DNA reported by Doležalová et al. (2002Doležalová et al. ( , 2003 and Sretenović Rajičić et al. (2008) in plant material obtained from germplasm collections. These authors detected notable differences in DNA content between the species and a highly unusual, puzzling genetic diversity in isozymes and AFLP profiles of L. livida that were attributed to misidentifications among germplasm accessions. Doležalová et al. (2002Doležalová et al. ( , 2003 pointed out the high genetic similarity of these materials to L. dregeana DC., a South African endemic related to L. sativa and L. serriola (Koopman et al. 2001). According to our results, most probably such erroneous species identifications are allowed by the striking uncertainty about L. livida delimitation.
unfortunately, the assumption that typical Lactuca virosa and L. livida morphotypes can readily cross has not been tested because of their strong ability to self-fertilize and the difficulty of emasculating the small and short-living florets (Mejías 1994). Nonetheless, bearing in mind these properties, it seems rather difficult to conceive the existence of independent tokogenetic systems or genotypic clusters within the complex and, thus, two differentiated biological species as the result of intrinsic constraints on genetic exchange or an effective isolation process in progress (Stace 1989;Mallet 1995;Wu 2001;Grant and Grant 2006;Koffi et al. 2010). The described biological uniformity does not suggest that we are dealing with cryptic species, i.e., two plants assemblages that are hardly distinguishable in terms of morphology but have different niches and diverging evolutionary trajectories. Molecular phylogenetic methodologies provide robust tools to reveal the occurrence of cryptic species (Bickford et al. 2007) and numerous cases in plants can be cited in this regard (e.g., Abdelaziz et al. 2011;Sokoloff et al. 2019, Li et al. 2020, Whittall et al. 2020. Commonly, the analysis involves the association of molecular data to ecological and/ or reproductive shifts with reduced morphological divergence, sometimes coupled to a geographical pattern. However, our study only provided a diffuse predominance of deeply lobed leaf shape in the western Iberian Peninsula without evidences of ecological differentiation or association to other traits. Despite this, a molecular phylogenetic study would be desirable to definitively and conclusively reject the hypothetical existence of distinct evolutionary lineages within L. virosa and to facilitate the correct identification of germplasm collections materials. Failure to detect any morphologically differentiated group also prevents the validation of any taxon at the infraspecific level within Lactuca virosa s.l. in the Iberian Peninsula (Velasco 1981). This argument is reinforced by our use of multivariate statistical methods, which have proven to be the most informative method of analysis at fine taxonomic levels, such as subspecies or variety (Michener and Sokal 1957;Sokal 1965;Stuessy 2009;Marhold 2011;Patten 2015). Moreover, the geographical and ecological differentiation that are major components in the recognition of infraspecific taxa (Stace 1989;Hamilton and Reichard 1992;Stuessy 2009) were virtually absent in our results. Among the characters tested, only leaf lobation shows a significant geographic pattern, but the distribution was only somewhat cohesive (Figure 4), which indicates that this variability is most likely due to general plasticity within the group. Thus, even a taxonomic treatment at the form rank is not recommended in this case (Davis and Heywood 1963;Stace 1989;Stuessy 2009).
The scarcity of evidence to support any taxonomic differentiation within Lactuca virosa in the Iberian Peninsula clearly indicates that the name Lactuca livida Boiss. & Reut. should be regarded as a synonym of L. virosa L. Given the results and discussion presented, we are somewhat surprised at the continued consideration of Lactuca livida Boiss. & Reut. as a valid taxon up to present times. It also seems striking that, despite the endemic distribution assigned (Feráková 1976(Feráková , 1977, it has never been included in Spanish threatened species lists (i.e., Bañares et al. 2004 and subsequent addenda). We wonder if botanists have long regarded L. livida as a taxon of doubtful validity, but none has been sufficiently concerned with the taxonomic problem to address it until now.

Taxonomic conclusions and implications for conservation
Our morphological analysis clearly shows that Lactuca virosa s.l. is not split into two phenetic species (Mallet 1995(Mallet , 2020Sites and Marshall 2004) on the Iberian Peninsula. This result, together with the biological knowledge of the group (Mejías 1993(Mejías , 1994, suggests no taxonomic differentiation processes under the biological, cladistic or phylogenetic species concepts. The absence of clear geographical and altitudinal differentiation also prevents the validation of any taxa at infraspecific rank, as some authors have proposed (Velasco 1981;Blanca and Cueto 2009). Thus, the assemblage of plants that, according to the literature (Boissier 1845;Feráková 1976Feráková , 1977Velasco 1981;Cueto 1985, 2009), could be identified as Lactuca livida Boiss. & Reut. or L. virosa subsp. livida (Boiss. & Reut.) Ladero & Velasco lacks any taxonomic identity, and these names must be strictly considered synonyms of L. virosa L. Although achieving taxonomic stability of classifications is a utopic goal because of the very nature of species (Agapow and Sluys 2005; Padial and de la Riva 2006), we hope that present considerations from diverse species concept perspectives, will contribute to a solid classification in the genus Lactuca. Moreover, we aim to facilitate Lactuca germplasm conservation, not only by avoiding unnecessary redundancy Doležalová et al. 2003Doležalová et al. , 2004 but also by preventing inconsistent identifications (Sretenović Rajičić et al. 2008).
The present work is a clear example of how morphometric analysis can be helpful to address troublesome taxa (Valcárcel and Vargas 2010;Jiménez-Mejías et al. 2014).
Here, morphometrics provided a powerful, straightforward tool for detecting false morphological diagnostic characters, thus enabling the synonymization and disambiguation of taxa (Pometti et al. 2007;Kougioumoutzis et al. 2016;Vogel Ely et al. 2018;Duan et al. 2019). Consequently, although poorly appreciated at present (Wiens 2007), morphometrics should be viewed as a useful approach to address problems of taxonomic instability and even inflation (Knapp et al. 2005), in particular when combined with ecological and biological information.