Spatial analysis
Occurrence records were obtained from 287 observations from the Global Biodiversity Information Facility (211 occurrences validated by morphology, http://www.gbif.org/) and data from herbarium specimens (76 occurrences). For modelling purposes, the dataset was reduced to keep only georeferenced data.
We used environmental niche modelling (ENM) approaches to reconstruct the potential distribution of the four main clades uncovered in the phylogenomic analyses (see Results) under current and past climatic conditions using the maximum entropy algorithm implemented in the R package ‘Maxent’ (Phillips, 2021). Nineteen bioclimatic variables were extracted from the Bioclim dataset, provided by WorldClim 1.4 in a GIS-based raster format (2.5 min resolution). The correlations between environmental variables were determined with a Pearson’s correlation matrix and subsequent realization of a dendrogram cluster for its visualization (Supplementary Data Figure S1). We selected a different set of uncorrelated variables for each geographical region with a high percentage of contribution (PC): bio4 (temperature seasonality), bio8 (mean temperature of wettest quarter), bio9 (mean temperature of driest quarter), bio15 (precipitation seasonality) and bio16 (precipitation of wettest quarter) for the circum-Mediterranean region; bio3 (isothermality), bio6 (min temperature of coldest month), bio8 (mean temperature of wettest quarter), bio14 (precipitation of driest month), bio15 (precipitation seasonality) and bio16 (precipitation of wettest quarter) for the Eastern Mediterranean region; bio3 (isothermality), bio4 (temperature seasonality), bio16 (precipitation of wettest quarter) and bio18 (precipitation of warmest quarter) for the Macaronesian region. ENM analyses were carried out under current climatic conditions, and projected to climatic conditions of the Mid Holocene (MH, about 6,000 years ago), the Last Glacial Maximum (LGM, about 22,000 years ago; Braconnot et al., 2007) and the Last Interglacial (LIG, ~120,000 - 140,000 years BP; Otto-bliesner et al., 2006) using the palaeoclimatic ‘Community Climate System Model’ (CCSM; Gent et al., 2011). Layers were cropped to represent the distribution range of each phylogenetic group (i.e., DC1, DC2 and DC3 for D. communis , and D. orientalis ; see Results) to maximize the reliability of the results and discard false occurrences using the package ‘raster’ (version 3.5-15; R v.4.0.5). We used the Schoener’s D and Hellinger’s I indices as implemented in ENMtools v.1.0.4 to evaluate niche overlap (Warren et al., 2008, 2010). Equivalence and similarity tests with 1,000 replicates were carried out to assess if the overlap between ENMs is higher than expected under randomized ENMs.
A multivariate ordination analysis (principal component analysis, PCA) was carried out using uncorrelated bioclimatic variables obtained from WorldClim using the packages ‘ade4’, ‘factoextra’, ‘magrittr’, ‘dismo’ and ‘HH’ in R v.4.0.5. The correlation analysis was performed with a Pearson correlation matrix and the subsequent visualization and selection of variables using a dendrogram cluster. The selection of uncorrelated variables with the highest contribution to the PCA were 10: bio1 (annual mean temperature), bio2 (mean diurnal range), bio3 (isothermality), bio7 (temperature annual range), bio8 (mean temperature of wettest quarter), bio9 (mean temperature of driest quarter), bio10 (mean temperature of warmest quarter), bio12 (annual precipitation), bio15 (precipitation seasonality) and bio19 (precipitation of coldest quarter).