INTRODUCTION

Accurate species delimitations are crucial for taxonomy and systematics, but also for any field that relies on biodiversity measures, such as macroecology and conservation (^{Agapow et al., 2004}). In this way, species provide data that are used for example to test hypotheses about large -scale ecological patterns, identify biodiversity hotspots, study patterns of threat and vulnerability, and rank areas for conservation (Isaac et al., 2004; Gippoliti & Amori, 2007). In order to achieve reliable delimitations for taxa, this poses the challenge for taxonomists to revise and test existing species hypotheses reflected in the names of Linnaean classification.

Xerodraba Skottsb. is a genus of Brassicaceae (Cruciferae) comprising five species [X. colobanthoides Skottsb., X. glebaria (Speg.) Skottsb., X. lycopodioides (Speg.) Skottsb., X. monantha (Gilg ex Kuntze) Skottsb., and X. patagonica (Speg.) Skottsb.] that inhabit southern Argentinean Patagonia, with two of them (X. lycopodioides and X. patagonica) also present in southern Chile (^{Salariato et al., 2015}a). The genus is principally differentiated from the remaining South American Brassicaceae by having dense cushions, solitary flowers terminating individual stems, imbricate scale-like leaves, and latiseptate silicles (Boelcke & Romanczuk, 1984; Al-Shehbaz, 2012). Xerodraba is currently included within Eudemeae Al-Shehbaz, Warwick, Mummenhoff & Koch, a tribe of 11 genera and 43 species distributed along the Andes from central Colombia at ~5° N to southern Argentinean Patagonia at ~55° S (Al-Shehbaz et al., 2023). It is phylogenetically related to the southern genera Onuris Phil. and Alshehbazia Salariato & Zuloaga (Salariato et al., 2015b, 2022), and molecular data indicate its early divergence in the Late Miocene-Pliocene (Salariato et al., 2022). Xerodraba originally included eight species (Skottsberg, 1916), though only the five listed above are currently recognized (Salariato et al., 2015a; Al-Shehbaz et al., 2023). Its presence in southern Patagonia is scattered and sparse, and the genus is poorly represented in the herbaria of the world, with some species known from only few specimens (e.g., X. glebaria) or even from just the type collection (e.g., X. monantha). Nowadays, X. patagonica, is the species with the largest geographic range in the genus, distributed in sympatry with the other species (with exception of X. glebaria), and differentiated mainly by having oblong to ovate and densely imbricate leaves (Salariato et al., 2015a; Al-Shehbaz et al., 2023). It includes two subspecies: X. patagonica (Speg.) Skottsb. subsp. patagonica and X. patagonica (Speg.) Skottsb. subsp. pycnophylloides (Speg.) Salariato & Al-Shehbaz. Originally, X. pycnophylloides (Speg.) Skottsb. was distinguished from X. patagonica by the length of its flowers and fruits (Skottsberg, 1916). However, subsequent authors pointed out a great similarity between the two species and were unable to distinguish them without their flowers (Boelcke & Romanczuk, 1984; Al-Shehbaz, 2012). Later, Salariato et al. (2015a) recognized that the differentiation between these two species is exclusively based on continuous variation mainly of the sepals, petals, and anther length. Following the phylogenetic species concept (PSC) of Nixon & Wheeler (1990) used in the modified population aggregation analysis (PAA) (Henderson, 2004, 2005), they recognized a single species with two subspecies. However, for an adequate taxonomic ranking of these names, more studies analyzing different types of data are still needed. Therefore, we aim here to shed light on the identity of these two names by contrasting results obtained from different phylogenetic, ecological niche, and morphological analyses.

MATERIALS AND METHODS

Phylogenetic analyses

For molecular analyses, we downloaded ITS (nrDNA dataset), trnL-F, trnH-psbA and rps16 intron (plastid DNA dataset) sequences of Xerodraba species from GenBank, mainly from ^{Salariato et al. (2015}b, 2022). We also generated 36 new DNA sequences for nine new specimens using Sanger sequencing (new sequences were submitted to GenBank, https:// www.ncbi.nlm. nih.gov/genbank/; voucher information and GenBank accession numbers are provided in Appendix S1 in the Supporting Information). For X. patagonica subsp. pycnophylloides, we were only able to include a single specimen from Salariato et al. (2015b), since herbarium material failed PCR amplification, and we were unable to find new populations in subsequent field trips. This fact prevents us from testing the reciprocal monophyly of both subspecies and, therefore, in the phylogenetic studies we focused only on the inclusion of X. patagonica subsp. pycnophylloides within X. patagonica subsp. patagonica (hereafter subsp. pycnophylloides and patagonica, respectively). This single specimen exhibits all morphological traits associated with subsp. pycnophylloides (see results of morphological analysis in this work - Fig. 3, and morphological analyses in Salariato et al., 2015a), and its collection locality is congruent with that of several specimens of subsp. pycnophylloides collected in southeastern Santa Cruz province (Salariato & Al-Shehbaz et al. 2015; Al-Shehbaz et al., 2023). Therefore, we consider for this study that this specimen constitutes a good representative of this subspecies. Protocols for extraction, amplification, and DNA sequencing follow Salariato et al. (2020). Sequences were assembled and edited using Chromas Pro v.2.1.10 (Technelysium Pty Ltd, Brisbane, Australia), which was also used to check for the presence of single peaks in the chromatograms, especially for ITS sequences. As outgroups, we include sequences of Onuris graminifolia Phil. and Alshehbazia hauthalii (Gilg & Muschl.) Salariato & Zuloaga. The final dataset included 19 accessions representing four Xerodraba species, while X. colobanthoides and X. monantha could not be sequenced. Alignments were first generated with MUSCLE v.3.8.31 (Edgar, 2004), and then manually checked and improved where necessary using Bioedit v.7.2.5 (Hall, 1999) . Aligned matrices and all other supplemental data are available from the Supporting Information and Figshare Digital Repository (doi: 10.6084/m9.figshare.23804469). Best-fit models of nucleotide evolution were identified using the Akaike information criterion corrected for small sample sizes (AICc) implemented in jModeltest2 v.2.1.6 (Darriba et al., 2012): TIM2ef+I (ITS), TIM2 (trnL-F), TPM2uf+I (trnH-psbA), F81 (rps16 intron).

First, we estimated molecular phylogenies of Xerodraba accessions to explore the phylogenetic placement of subsp. patagonica and subsp. pycnophylloides. Phylogenies were calculated for the ITS, cpDNA, and combined ITS+cpDNA datasets using Bayesian inference in MrBayes v.3.2.7a (Ronquist et al., 2012), setting the number of substitution types to “mixed” (which results in the Markov chain sampling over the space of all possible reversible substitution models), and rates according to the recovered best- fit models. Two simultaneous analyses, starting from different random trees and with four Markov Monte Carlo chains, were run for 10 million generations, sampling every 10,000 generations to ensure the independence of the successive samples. The convergence and effective-sample size were checked with the average standard deviation of split frequencies (ASDSF) < 0.01, the potential scale reduction factor (PSRF) ~ 1, and the effective sample size (ESS) for all parameters - 200. After checking for convergence and effective-sample size, the first 250 trees (25 % of total trees) were discarded as burn-in, and the remaining samples of each run were combined and used to calculate the 50 % majority rule consensus tree. Incongruences between ITS and cpDNA data were also explored in a filtered supernetwork calculated with SplitsTree 4.14.8 (Huson & Bryant, 2006) using the Z-closure algorithm (Huson et al., 2004), 1502 Bayesian posterior trees of each nuclear and plastid dataset, and filtering the splits to show only those present in a minimum of 30 % input trees.

Second, we employed the Bayes factor-delimitation approach (BFD) (Grummer et al., 2014) to test the support for the hypothesis that both subspecies represent the same species (vs. the hypothesis that they represent two different species). Different hypotheses were tested with the ITS and cpDNA data using the multispecies coalescent model implemented in *BEAST extension (Heled & Drummond, 2010) of BEAST v1.8.4 (Drummond et al., 2012). For the single-species hypothesis, accessions of both subsp. patagonica and pycnophylloides were included under the same species, while for the alternative hypothesis, specimens corresponding to each subspecies were placed under different species. Analyses in *BEAST were conducted with the nucleotide substitution models selected in jModeltest2, an uncorrelated lognormal clock model (UCLN), a Yule process for the species tree prior, and the piecewise linear with constant root for the population size model assigned to each locus. For each hypothesis, three runs were conducted using 50 million generations and sampling every 12,500.

The BF was calculated using MLE obtained both by path-sampling (PS; Lartillot & Philippe, 2006) and stepping-stone sampling methods (Xie et al., 2011) with 100 steps of one million generations each and α = 0.3. The 2log(BF) = 2[log(ML1)-log(ML2)] was calculated to compare the competing hypotheses following criteria of Kass & Raftery (1995): 2log(BF) = 0-1 “not worth more than a bare mention”, 2log(BF) = 1-2.6 “positive” support, 2log(BF) = 2.6-4.4 “strong” support, and 2log(BF) > 4.4 “decisive” support for model 1 over 2.

Climatic Niche analyses and geographic range overlap

To study and compare climatic niches of patagonica and pycnophylloides subspecies, we used occurrence data for Xerodraba obtained from the examination of specimens deposited in different herbaria and field collections (specimen vouchers in Table S1, Supplementary Material). All records were previously mapped using QGIS v3.26.3 “Buenos Aires” (Quantum GIS Development Team, 2022) for visual inspection, and in cases of specimens with no GPS coordinates but exact locality names, records were georeferenced using Google Earth Pro v7.3.6.9345 (https://www.google.com/intl/en/earth/). After removing duplicates and occurrences closer to 30 arc-seconds (~ 1 km), we obtained a total of 77 records for Xerodraba, of which 23 corresponded to subsp. patagonica, 11 to subsp. pycnophylloides, 5 to X. colobanthoides, 3 to X. glebaria, and 34 for X lycopodioides. These records corresponding to the complete genus were used to delimit the study area (see below). Data on the current climatic conditions were extracted from the CHELSA 2.1 climatic dataset (Karger et al., 2017, 2021). Values for 19 bioclimatic variables were extracted at a resolution of 30 arc-seconds from the study area, defined by a minimum convex polygon enclosing all specimen records for the genus with a 100-km buffer zone (Fig. S1 from Supplementary Material). Data extraction and manipulation were done using the R (R Core Team, 2022) packages adehabitatHR v.0.4.20 (^{Calenge, 2006}), raster v.3.6.20 (Hijmans, 2023), sp v.1.6.0 (Pebesma & Bivand, 2005; Bivand et al., 2013) and maptools v.1.1.6 (Bivand & Lewin-Koh, 2022).

Dimensionality of environmental space was then reduced by performing a PCA-env approach (Broennimann et al., 2012) with the R package ade4 v. 1.7.22 (Dray & Dufour, 2007), in which a principal component analysis (based on a correlation matrix) is calibrated on the entire environmental space within the study area. For subsequent analysis, we retained the first five PCs (principal components), which accounted for 97.75 % of total climatic variation (Table S2, Supplementary Material).

For niche comparison between both subspecies we applied the approach of Broennimann et al. (2012), in which environmental space was gridded in 100 x 100 cells (each cell corresponds to a unique vector of the available environmental conditions in the study area) and a kernel density function was applied for smoothing the density of occurrences for each of the cells in the environmental space (see Broennimann et al., 2012; Di Cola et al., 2017 for details on the methodology and the kernel density estimator). Density grids for each subspecies were used subsequently to compute the niche overlap by means of the Schoener’s D statistic (Schoener, 1968) and the niche expansion (representing the amount of niche of subsp. pycnophyllodes absent in niche of subsp. patagonica), stability (amount of niche of subsp. pycnophyllodes included in niche of subsp. patagonica), and unfilling (amount of niche of subsp. patagonica absent in niche of subsp. pycnophylloides) indices (Di Cola et al., 2017). Although these indices are frequently used to compare the niche of the same species distributed in different regions or at different times (e.g., in its native and non-native range, Guisan et al., 2014), in this case they were calculated to explore the hypothesis that climatic niche of subsp. pycnophylloides is included within that of subsp. patagonica, as they represent the same species, hence the order of the comparisons made. Significance of these indices was assessed using the similarity test (Warren et al., 2008) and 9999 replications with a two-tailed test, in where overlap = higher, expansion = lower, stability = higher, unfilling = lower support niche similarity hypothesis; and overlap = lower, expansion = higher, stability = lower, unfilling = higher supports niche divergence hypothesis.

All niche comparisons in environmental space were performed using the ecospat package v3.5 (Broennimann et al., 2023).

Alternatively, we estimated and compared the distribution ranges of both subspecies using species distribution modeling. We use the same first five PCs of the environmental space for the study area at a resolution of 30 arc sec as predictor variables (which captured 97.75% of the total climatic variation contained in the 19 original bioclimatic variables for the study area), applying the ensemble of small modes (ESM) technique, which is suitable when the number of occurrences is small (Breiner et al., 2015). For this purpose, these five covariates were used in bivariate models (only two predictors at a time) considering all possible combinations (10 covariate combinations) and calibrated by using MaxEnt (Elith et al., 2011), since this technique reports good performance for ESM (Breiner et al., 2018). The total number of pseudo-absences/background points sampled for each species was set to 10000, weighting presences equally to pseudo-absences (prevalence = 0.5). For each bivariate model, we conducted ten cross-validations (80 % - 20 %), and predictive performance was evaluated using the true skill statistic (TSS) and the area under the receiver-operating characteristic curve (AUC) metrics. For each species we obtained 100 predictive models (10 models × 10 cross-validations), which were used to generate an ensemble model based on the weighted sum of probabilities of each model as evaluated by the TSS scores, dropping models with TSS < 0.5. Then, for each subspecies, we calculate variable contribution of each PC as the ratio between sum of weights of bivariate models (based on AUC scores) where a focal variable was used and sum of weights of bivariate models where the focal variable was not used (corrected for the number of models with or without the focal variable) . In this way, we compared variables (climatic PCs) that contributed most to the modeling of subsp. patagonica vs. that of subsp. pycnophylloides. In the last step, to convert continuous ensemble predictions of each subspecies into binary (presence/absence) maps, we use the maximum training sensitivity plus specificity threshold (MaxTss, Liu et al., 2013).

The distribution range of subsp. patagonica and subsp. pycnophylloides was compared by counting the number of shared and unique cells for the presence of each subspecies. Additionally, modeling was also carried out, including all the records of both subspecies as a single species, so that we can compare the distribution range obtained for X. patagonica (subsp. Patagonica + subsp. pycnophylloides) with that of subsp. patagonica. Species distribution modeling and geographic-range comparison were carried out with the ecospat and biomod2 v3.5.1 (Thuiller et al., 2021) packages.

Morphological analyses

We first performed a preliminary analysis of morphological variation based on ^{Salariato et al. (2015}a) over 34 specimens of Xerodraba patagonica (24 for subsp. patagonica and 10 for subsp. pycnophylloides). The preliminary traits studied were leaf length, maximum width of leaf, trichome length on leaf margins, sepal length, petal length, maximum width of petal, anther length of median stamens, style length in mature fruits, length of fruit valve, and maximum width of fruit valve. For subsequent analyses, we selected four quantitative characters that showed variation between the two subspecies: leaf length (mm), sepal length (mm), petal length (mm), and anther length of median stamens (mm), all of which with p < 0.01 in the Wilcoxon rank-sum test when comparing values for subsp. patagonica vs. subsp. pycnophylloides (Fig. S2 from Supplementary Material and Table S1). Then, for these four traits, we conducted a PCA (Principal Component Analysis) with the R package ade4 v. 1.7.22 (Dray & Dufour, 2007) based on the correlation matrix, obtaining 95.69 % of the total variation in the first three components. These first three PCs, for which sepal and petal length dominated in PC1, leaf length in PC2, and anther length in PC3 (Table S3, Supplementary Material), were used in posterior analyses. The efficiency of each PC in distinguishing these two subspecies was evaluated using the random forest (RF) method (Breiman, 2001) with the packages randomForest v4.7.1.1 (Liaw & Wiener, 2002) and caret v6.0-93 (Kuhn, 2022), with a maximum of 500 trees and evaluating the performance of the model using leave-one-out cross-validation and the Cohen’s Kappa statistic. We then studied the relation between geographic distribution, morphological variation, and classification of both subspecies, for which: (1) we first evaluated by RF the performance of a classification model based only on geographic distribution (i.e., subsp. ~ lat + long), (2) generated independent predictions for PC1 using leave- one-out cross-validations from a second-order trend surface linear model (i.e., PC1 ~ lat + long + lat^2 + long^2 + lat:long interaction), which were later used in RF to classify subspecies, evaluating performance in the same way as with the original PCs. Following a similar framework, relation between climate niche, morphological variation, and classification of both subspecies were analyzed, (1) evaluating by RF the performance of a classification model based only on the five climatic components of the environmental space (i.e., subsp. ~ climatic PC1 + climatic PC2 + climatic PC3 + climatic PC4 + climatic PC5), and (2) indirectly generating predictions of morphological PC1 with a linear model including linear and quadratic terms of climatic PCs as predictors, and then using the predicted morphological PC1 in RF to classify subspecies, evaluating performance with leave -one-out cross validation and the Cohen’s Kappa statistic.

Threat assessments

Preliminary threat assessments for recognized taxa were based on the IUCN Red List categories and criteria v3.1 (IUCN, 2012) following the IUCN guidelines v14 (IUCN, 2022). Conservation status were determined under criterion B, which is based on geographic ranges (IUCN, 2012; 2022) and suitable for estimating conservation status when the distribution of taxa is only known from georeferenced herbarium specimens with limited information on abundance and potential continuing decline (Nic Lughadha et al., 2018). For categorization we used extent of occurrence (EOO) (subcriterion B1; IUCN, 2022) calculated with a minimum convex polygon around occurrence points. For condition “b” we evaluated vulnerability and potential deterioration of the habitat where taxa grow. All IUCN parameters needed for assessments of taxa under criterion B were calculated using the R package ConR v1.3.0 (Dauby, 2020).

RESULTS

Phylogenetic analyses

Phylogenies recovered by the different analyses show Xerodraba composed of two main lineages, of which one includes X. lycopodioides, and the other includes X. patagonica (both subspecies) and X. glebaria (Fig. 1). In these phylogenies subsp. pycnophylloides was included within the X. patagonica lineage. Likewise, filtered supernetwork estimated with the ITS and cpDNA trees recovered specimen grouping similar to that of the phylogenetic trees (Fig. 1). When species trees and Bayes factor values corresponding to the different hypotheses (single-species vs. two-species hypotheses) were estimated, we found that neither hypothesis prevailed over the other (Table 1).

Table 1 Bayes factor species delimitation results using the multispecies Coalescent model in *BEAST v1.8.4. First model “PAT, PYC” corresponds to the two-species scenario (X. patagonica and X. pycnophylloides) while the second model “(PAT PYC)” represents the single species scenario (X. patagonica subsp. patagonica + X. patagonica subsp. pycnophylloides). Bayes factors are calculated between the two-species model vs. the single species model. Values correspond to the mean of three replicates. Marginal likelihood was estimated both by path sampling (PS) and stepping stone (SS) approaches. 

Fig. 1 Phylogenetic placement of X. patagonica subspecies within Xerodraba. A-C, Bayesian 50 % majority-rule consensus trees from 1502 trees generated by Bayesian inference with MrBayes v3.2.7a showing the phylogenetic placement of X. patagonica subsp. patagonica (identified as patagonica) and of X. patagonica subsp. pycnophylloides (identified as pycnophylloides). A, nrITS dataset. B, plastid DNA dataset (trnL-F/trnH-psbA/rps16 intron). C, combined nrITS+cpDNA. Values on branches correspond to Bayesian posterior probability (%). Units of branch length are proportional to nucleotide substitutions per site. D, filtered supernetwork generated in SplitsTree v4.14.8 using 1502 Bayesian posterior trees per each nuclear ITS and plastid DNA dataset, and filtering the splits to show only those present in a minimum of 30 % input trees). In all cases the asterisk indicates position of subsp. pycnophylloides. Color version at https://www.ojs.darwin.edu.ar/index.php/darwiniana/article/view/1149/1317 

Fig. 2 Niche and range comparisons between subspecies of X. patagonica. A-C, niche comparisons along the environmental (E) space produced by the two main axes of the PCA-env. For each subspecies, the gray to black shading represents the grid-cell density of occurrence (black being the highest density). The dashed line represents 50 % of the available environment and the solid line represents 100 %. A, X. patagonica subsp. patagonica. B, X. patagonica subsp. pycnophylloides. C, niches for both subspecies, where the (E)-space colored with green corresponds to the portion of the subsp. patagonica niche not occupied by subsp. pycnophylloides, while the space colored with blue shows the niche space shared by both subspecies. D-G, distribution modeling for both subspecies using ESMs and MaxEnt. D-E, predicted suitable climatic conditions (logistic output) from the MaxEnt model for subsp. patagonica (D) and subsp. pycnophylloides (E) using the first five principal components of the PCA-env as climatic variables. Dots on the maps correspond to the occurrences of each subspecies. F, binary (presence/absence) distributions maps for both subspecies derived from the ESMs outputs using the maximum training sensitivity plus specificity as threshold. Yellow cells indicate presence only of subsp. patagonica, red cells presence only of subsp. pycnophylloides, and green cells indicate presence of both subspecies. G, binary map obtained when occurrences of the two subspecies were included together as a single species. Color version at https://www.ojs.darwin.edu.ar/index.php/darwiniana/article/view/1149/1317 

Climatic Niche analyses and geographic range overlap Niche comparisons in the environmental space showed that overall, the niche of subsp.

pycnophylloides is included within that of subsp. patagonica (Fig. 2A-C). In particular, the result was reflected in values of the expansion (amount of niche of subsp. pycnophyllodes absent in niche of subsp. patagonica) and stability (amount of niche of subsp. pycnophyllodes included in niche of subsp. patagonica) metrics, which were close to 0 and 1, respectively (Table 2). Moreover, when the similarity test was significant (p ≤ 0.05), it indicated the presence of niche conservatism (overlap = higher, expansion = lower, stability = higher, unfilling = lower), but no significance was obtained for niche divergence with any of the climatic PCs (Table 2). On the other hand, species distribution models obtained by ESM for the two subspecies (subsp. patagonica: TSS = 0.629, AUC = 0.834; subsp. pycnophylloides: TSS = 0.634, AUC = 0.883) agreed in recovering

Table 2 Values for the different metrics obtained from the comparisons between climatic niches of subsp. patagonica and subsp. pycnophylloides in the environmental space. Numbers correspond to the observed values, while in brackets is included the p corresponding to 9999 replicates of the two-tailed similarity test. In all cases, significance (in bold) was associated with the niche conservatism hypothesis (i.e., overlap=higher, expansion=lower, stability=higher, unfilling=lower). 

PC3 and PC4 with higher contribution, and PC1 and PC2 with lower contribution, while models of both subspecies differed in the contribution of PC5 (higher in subsp. patagonica and lower in subsp. pycnophylloides) (Table 3). Presence-absence maps obtained with ensemble models and the MaxTss threshold showed that range of subsp. pycnophylloides is largely overlapping with those of subsp. patagonica, with 99.4 % of the range of the former included in the range of the latter (Fig. 2D-E). When comparing distribution area of subsp. patagonica with the distribution area of X. patagonica (for the latter case modeling included the records of both subspecies: TSS = 0.592, AUC = 0.841), we recovered that 85.3% of the species range (subsp. pycnophylloides + subsp. patagonica) was included in the subspecies range (subsp. patagonica) (Fig. 2F).

Morphological analyses

From the three morphological PCs analyzed (Fig. 3), only PC1, mainly related to the length of sepals and petals (Fig. 3, Table S3), showed high performance to distinguish between both subspecies (Cohen’s Kappa statistic = 1, Table 4). When geographic position was used as predictor, low kappa values were obtained, both when it was included directly as well as when its effect was included through PC1 (Fig. 3, Table 3). Similarly, lower kappa values were also recovered when the climatic PCs were used as predictors (either directly or indirectly through morphological PC1), showing the limitations of both types of data (geographic and ecological) in identifying these subspecies.

DISCUSSION

From the results of different types of data, we can highlight five main points: (1) molecular data do not support the separation of subsp. pycnophylloides as a distinct species from X. patagonica; (2) petal and sepal length effectively distinguishes the two subspecies, representing different portions of the continuous variation in the size of floral parts (larger in subsp. patagonica and smaller in subsp. pycnophylloides); (3) climatic niche of subsp. pycnophylloides does not differ from that of subsp. patagonica in any of its dimensions, by contrast, for some of these components significant niche conservatism is recovered; (4) both taxa show pronounced sympatry in their distributions along the Patagonian steppe; (5) morphological differentiation of both subspecies associated with the length of their floral parts does not seem to be directly related to either their geographic distribution or their ecological niche.

Fig. 3 Morphological variation in X. patagonica subspecies. A, biplot of the first two principal components extracted in the principal component analyses (PCA) using four characters: leaf length, sepal length, petal length, and anther length. Inertia ellipses correspond to one standard deviation (~ 68 % of variation) for each subspecies. Factor loadings are shown in Table S3 (Supplementary Material). B, variation of PC1 along the geographical space, where sizes of circles are proportional to their PC1 values. Numbers in A and B indicate position of specimens included in the phylogenetic analysis. Color version at https://www.ojs.darwin.edu.ar/index.php/darwiniana/article/view/1149/1317 

Table 3 Variable contribution of each variable (climatic PCs) for species distribution models of X. patagonica subspecies calculated as the ratio between sum of weights of bivariate models where a focal variable was used and sum of weights of bivariate models where the focal variable was not used. In bold the climatic PC with the highest contribution. 

Table 4 Cohen’s Kappa statistic values for the different models evaluated using the random forest technique and leave-one-out cross-validation. PC1m, PC2m, PC3m refers to the different morphological PCs. 

Morphological and molecular evidence (i.e., the high morphological similarity between the two subspecies (^{Salariato et al., 2015}a) and the phylogenetic placement of subsp. pycnophylloides within subsp. patagonica) suggest that populations/metapopulations associated with these two names belong to the same independent evolutionary lineage, also interpreted as the same evolutionary significant unit (ESUs; Moritz, 1994). Following the General Lineage Concept (de Queiroz, 1998, 1999, 2007), which treats species as separately evolving metapopulation lineages through time (Simpson, 1951; de Queiroz, 2007), these names should be included under the same species, focusing the question now on the validity of both subspecies. Beyond the diversity of species concepts, subspecies and varieties are usually defined as metapopulations with slight morphological or phylogenetic differences due to geographical or ecological differentiation (^{Reydon & Kunz, 2021}). This definition highlights two perspectives for delimiting infraspecific taxa:

(1) geographic subspecies, isolated by distance or landscape barriers; (2) ecological subspecies, occupying different ecological niches (Meikle, 1957). Results from climatic niche analyses and geographic range estimations show sympatry and niche similarity for both subspecies, therefore, we consider that there is no evidence to justify their maintenance as subspecies. In addition to the study of phylogenetic structure, morphometric and ecological niche analyses have been successfully used to study the delimitation of taxa in Andean and Patagonian species complexes (e.g., ^{Salariato et al., 2012}; Nicola et al., 2014; Ferrero et al., 2020). Here, continuous variation of petal, sepal, and anther length in subsp. pycnophylloides and subsp. patagonica appears to represent only extremes of the variation range of X. patagonica, and since it does not appear to be related to differentiation in their geographic distributions or ecological niches, it lacks support for the maintenance of both subspecies. Consequently, and based on this evidence, we have decided to unite both subspecies.

Regarding X. glebaria, molecular phylogenies place it within the X. patagonica lineage. However, X. glebaria is easily distinguished from the remaining species by its small, ovate leaves with margins fimbriate-ciliate all around (including the apex) with densely packed, simple trichomes up to 0.6 mm long, and by sepals ciliate at margins (^{Salariato et al., 2015}a; Al-Shehbaz et al., 2023). This species, so far represented by less than 10 specimens, and of which no new collections have been reported in the last 70 years, has been collected mainly in the south of Chubut province (Depts. Río Senguer and Sarmiento), in localities where X. patagonica is absent or at least it has not been found to date. Given the consistency and stability of morphological characters that distinguish X. glebaria, its non-overlapping distribution range with respect to that of X. patagonica, and the lack of new collections that provide more evidence of morphological and geographic variation, we prefer to maintain the species and not include it within the synonymy of X. patagonica until new evidence supports this decision. For analysis and discussion of the remaining species and synapomorphies of the genus, see Salariato et al. (2015a).

TAXONOMIC TREATMENT

Xerodraba patagonica (Speg.) Skottsb., Kongl. Svenska Vetensk. Acad. Handl., n.s. 56(5): 362, tab. 22, fig. 13. 1916. Braya patagonica Speg., Anales Soc. Ci. Argent. 47: 171. 1899. Eudema patagonica (Speg.) Gilg & Muschler, Bot. Jahrb. Syst. 42(5): 471. 1909. Type: Argentina. [Province of Santa Cruz]. Depto. Lago Argentino, Lago Argentino, Karr-aik, March 1898, C. Ameghino s.n. (holotype LPS-10432! in LP!; isotypes BAA! fragm. ex LP, UPS!).

Braya pycnophylloides Speg., Anales Soc. Ci. Argent. 47: 172. 1899. Eudema pycnophylloides (Speg.) Gilg & Muschler, Bot. Jahrb. Syst. 42(5): 471. 1909. Xerodraba pycnophylloides (Speg.) Skottsb., Kongl. Svenska Vetensk. Acad. Handl., n.s. 56(5): 362, tab. 22, fig. 15. 1916. Xerodraba patagonica (Speg.) Skottsb. subsp. pycnophylloides (Speg.) ^{Salariato & Al-Shehbaz, Phytotaxa 207(1): 61. 2015}, syn. nov. Type: Argentina. [Province of Santa Cruz]. Depto. Lago Argentino, Lago Argentino, hab in montains Karr-aik, prope Lago Argentino, March 1898, C. Ameghino s.n. (holotype LPS-19006! in LP, isotypes, B!, BAA!, UPS!).

Eudema microphylla Gilg & Muschler, Bot. Jahrb. Syst. 42(5): 472. 1909. Xerodraba microphylla (Gilg & Muschl.) Skottsb., Kongl. Svenska Vetensk. Acad. Handl., n.s. 56(5): 362, tab. 22, fig. 11. 1916. Xerodraba pycnophylloides var. microphylla (Gilg & Muschl.) O.E. Schulz (1924: 250). Type: Argentina. [Province of Santa Cruz] Depto. Güer Aike, Patagonia Austral, Río Gallegos, O. Nordensjöld A-56 (holotype B!; isotypes BAA! fragm. ex B!).

Plants forming compact cushions, ultimate branches 1-5 cm tall; caudex up to 2 cm in diam., branches covered with persistent leaves of previous years. Leaves densely imbricate, straight, oblong-ovate, firm, (1.6-)2 -3.5(-4) × 0.6-1.6 mm, glabrous except for ciliate margin, abaxially not carinate, midvein obscure, adaxially concave to nearly flat, not expanded at base, margin ciliate along proximal half with separated simple trichomes 0.05-0.1(-0.15) mm long, apex obtuse. Fruiting pedicels slender, (1.5-)2- 4 mm long, glabrous. Sepals oblong, 2.5-5.2 × 0.9-1.8 mm, slightly fleshy at center, erect-ascending, subapically with simple trichomes to 0.3 mm long, margin membranous or rarely minutely ciliate; petals white to pale yellow, narrowly spatulate, 4-10.8 × 1-1.7 mm, apex rounded; claw slender, unappendaged; filaments 2-5 mm long, not expanded at base; anthers oblong, 0.7-1.5 mm long; nectar glands tooth-like; ovules 4-6 per ovary. Fruits ellipsoid, (3.5 -)4-6 × 2-3.6 mm; valves leathery, carinate; style 0.3-0.6 mm long. Seeds 1.5-3 × 1-2 mm (Fig. 4).

Fig. 4 Xerodraba patagonica. A, habit. B, plant with flowers. C, plant with fruits. A-B from Zuloaga et al 14743 (SI), C from Zuloaga et al. 14104 (SI). Photographs by Diego L. Salariato. Color version at https://www.ojs.darwin. edu.ar/index.php/darwiniana/article/view/1149/1317 

Distribution and habitat. Xerodraba patagonica grows along Andean slopes and Patagonian steppe in southern Argentina (provinces of Chubut, Neuquén, and Santa Cruz) and Chile (Region of Magallanes and Antartica Chilena) at 20-1050 m a.s.l. Its presence in the Neuquén Province is based only on H.F. Comber 335 (E, K), assigned to Cerro Bayo, but with the province not assigned. Since there are other localities called “Cerro Bayo” in Santa Cruz where the species is distributed (e.g., in Dept. Lago Buenos Aires and Dept. Deseado), new collections will be necessary to confirm the presence of this species (and the genus) in Neuquén.

Phenology. This species flowers from October into March.

Chromosome number. 2n = 22 (for X. pycnophylloides) (Manton, 1932).

Conservation status. Preliminary threat assessment for X. patagonica based on the IUCN Red List categories and criteria v3.1 (IUCN 2012, 2022) under subcriterion B1(a,b) gives it the LC status. However, populations of this species, as well as those of others distributed along the Patagonian steppe, should be monitored, since this region suffers from intense livestock grazing, as well as the development of mining and hydrocarbon activities, which has accelerated desertification and soil-erosion processes (Mazzonia & Vazquez, 2009; ^{Chartier et al., 2013}; Funk et al., 2018).

Specimens examined

Argentina. Province not assigned: Cerro Bayo, 25 February 1925, Comber 335 (E, K). Chubut: Depto. Escalante, Pampa del Castillo, 13-X-1929, Donat 181 (BM, CAS, F, GH, K, MO, S, SI, UC, Z); Donat 182 (BAA); Pampa del Castillo, 16-XI-1975, Garrido 734 (BAA, BAB); Pampa del Castillo, 2-II-1932, Castellanos s.n. (SI-167213, BA-6011); Rada Tilly, golfo San Jorge, without collector (LPS-21541 in LP); Ruta Provincial 27 hacia Río Chico a unos 10 km del cruce con la Ruta Nacional 3, 28-X-2005, Cocucci et al. 3587 (SI). Santa Cruz: without locality, Burmeister 169 (BAB, SI); Depto. Corpen Aike, Ruta Nacional 3, Parque Nacional Monte León, 26-XI-2012, Zuloaga et al. 14111 (SI); Parque Nacional Monte León, sobre Ruta Nacional 3, 6-XI-2005, Cocucci & Sérsic 3705 (CORD); a 111,6 km desde el cruce de la RP 9 con RN 3, rumbo a Calafate, 23-I-2013, Barboza et al. 3712 (CORD); Depto. Deseado, Barda Las lajas, Palau Mt., Comber 260 (E, K); Depto. Güer Aike, 500 m de Ruta Nacional 3, entrada a Yacimiento Estancia La Maggie, 26-XI-2012, Zuloaga et al. 14104 (SI); Estancia las Vizcachas, ca. 2 km O de laguna Las Vizcachas, 23-I-1977, TBPA 2372 (BAA, HIP); Estancia Le Marchand sobre Ruta Nacional 3, 6-XI-2005, Cocucci et al. 3703 (SI); Güer Aike, sobre Ruta Nacional 5 a Calafate, unos 3 km de Güer Aike, 5-XI-2005, Cocucci et al. 3701 (SI); Hotel Esperanza, 16-III-1914, Iter Patagonicum 360 (SI); orillas del Río Coyle, Bardas Blancas, 30-IX-1950, Sleumer 919 (CONC, US); por Ruta 7, a 2.5 km del cruce con Ruta Nacional 5, en dirección a Estancia Camasu Aike, 1-II-1995, Fortunato & Elechosa 4919 (BAB); Río Gallegos, Cerro Tres Hermanos, 2-XII- 1950, Sleumer 930 (CONC, US); Ruta Nacional 40, 15 km N de la Esperanza hacia Lago Argentino, 27-XI-1963, Correa 2868 (BAA); Estancia Güer Aike, 28-XI- 1950, Sleumer 891a (BAA, G, US); Río Gallegos, Kalela 1888 (H); 1897, Nordenskjöld s.n. (UPS); 10-XII-1975, TBPA 397 (HIP, SI); Río Gallegos, Estancia Los Pozos, 5-XII-1975, TBPA 129 (BAA); Río Gallegos, laguna La Leona, 7-XII-1975, TBPA 266a (BAA, BAB, SI); Depto. Lago Argentino, Cerro El Calafate, en la cumbre, 20-XII-1950, Sleumer 1201b (BAA); Cerro Huiliches, Laberinto 2, Olla Chingue, 15-XII-2001, Guerrido & Fernández 446a (SI); Cerro Huiliches, 23-XI-2012, Zuloaga et al. 14035 (SI); Zuloaga et al. 14043 (SI); 19-XI-2013, Zuloaga et al. 14734 (SI); Lago Argentino, Boulders Cerro Calafate, 25-XI-2007, Guerrido et al. 729 (SI); El Calafate, Cerro Frías, 19-XI-2013, Zuloaga et al. 14743 (SI); camino de Mata Amarilla a Estancia Australasia, 10 km al sur de la Estancia Bajada de la Oriental, 16-III-1977, León 2189 (BAA); Depto. Lago Buenos Aires, Patagonia prope Lago Buenos Aires, 15-XII-1902, Koslowsky s.n. (BAA); Depto. Río Chico, La Bajada, 1-III-1914, Iter Patagonicum 377 (SI). Chile. XII Magallanes and Antartica Chilena: Prov. Magallanes, Punta Arenas, Península de Brunswick, Reserva Forestal Laguna el Parrillar al norte, 18-II-1993, Henríquez & Pisano 136 (CONC, HIP); Península de Brunswick, Laguna el Parrillar al norte, 14-I-1982, Dollenz 947 (HIP); Cumbre del cerro próximo al monte Fen, 19-II-2003, Domínguez 1 (CONC).