We examined the realized, retrospective and predicted niches separately for J. phoenicea s.s., J. turbinata and J. canariensis and for the entire J. phoenicea complex. Additionally, we analysed separately the data from the groups of localities of J. turbinata detected in the genetic study. Juniperus phoenicea complex (ENT), J. phoenicea s.s. (PHO), J. canariensis (CAN) and J. turbinata (TUR); four different parts of J. turbinata’s geographic range: Atlantic (TURAT), central and west Mediterranean (TURCM), east Mediterranean (TUREM) and the Arabian Peninsula (TURAR).

The model was evaluated with ENMeval R software. We used 10 k-fold spatial partitions for each species presence record and evaluated models with the following feature classes: linear, quadratic, hinge, product and threshold, and the following values of regularization multipliers: 0.5, 1, 1.5, 2, 2.5, 3, 3.5, and 4.

Bioclimatic variables were obtained from the WorldClim 2.1 database (http://worldclim.org/) and Paleoclim (http://www.paleoclim.org/). We used: Last Interglacial period (LIG), Last Glacial Maximum climate (LGM) and Mid-Holocene (MH), present and the predicted warming in scenarios RCP2.6 and RCP8.5.

Data on the occurrence of the J. phoenicea complex was obtained from the Global Biodiversity Information Facility (GBIF.org), the literature, herbaria and the authors’ field notes.