Centaurea Subsect. Phalolepis (Compositae, Cardueae): A Case Study of Mountain-Driven Allopatric Speciation in the Mediterranean Peninsulas

Centaurea subsection Phalolepis has been thoroughly analyzed in previous studies using microsatellites in four centers of speciation: Anatolia, Greece, the Italian Peninsula and the Iberian Peninsula. Evidence suggests a correlation between taxon diversity and mountains. This group constituted a good case study for examining the mountain–geobiodiversity hypothesis (MGH), which explains the possible reasons for the many radiations occurring in mountains across the world. We combined all the datasets and carried out analyses of their genetic structure to confirm the species of subsect. Phalolepis are grouped according to a geographic pattern. We then checked whether climatic fluctuations favored the “species pump” hypothesis in the mountains by using the Climatic Stability Index (CSI). Finally, the relief of the terrain was tested against the rate of allopatric speciation by region by means of Terrain Ruggedness Index and environmental gradients through our new Climate Niche Breadth Index. Our results supported the MGH hypothesis and confirmed that the main triggers, namely altitudinal zonation, climatic oscillations and rugged terrain, must be present for the development of a radiation.


Introduction
Mountains are often linked to high levels of biodiversity [1,2]. In recent times, a study focused on the Tibeto-Himalayan region (THR) led to the definition of the "mountaingeobiodiversity hypothesis" or MGH to explain its species richness [3]. These authors proposed that three conditions are required to maximize the impact of mountain formations and surface uplift on regional biodiversity patterns, which are key to the origin of montane biodiversity. These are (i) complete altitudinal zoning with lowland, montane and alpine zones; (ii) climatic fluctuations that enable mountains to act as "species pumps"; and (iii) high-relief terrain with large environmental gradients that favor both species persistence in suitable refugia and allopatric speciation driven by geographic barriers. This hypothesis is very attractive and fits partly our own results in the Tibet-Himalaya-Hengduan region for the genus Saussurea (Herrando et al., pers.comm.). However, extrapolation of this hypothesis to other mountain systems remains unclear, despite the claims about this extrapolation made by [4].
In the Mediterranean environment, we do not find extremely rugged mountainous regions, but there might be some with the conditions proposed by [3] needed for high levels species is the highest, and the studies confirmed that they are genetically well-defined species [10,11]. The number of species decreases in an east-west gradient as the roughness of the mountains diminishes: in the Italian Peninsula, genetic studies do not support the existence of the microspecies described in the south of the peninsula [12]. As for the Iberian Peninsula, our study ruled out the presence of allopatric speciation and suggested that all speciation events in the group occurred via hybridization [13].
With this background, the most consistent hypothesis is that the southern Mediterranean mountains have promoted allopatric speciation in Phalolepis, favoring population isolation and allopatry, likely following the MGH model [3]. To test this hypothesis, we pooled the genetic data obtained for the same set of microsatellites from the four Phalolepis speciation centers and conducted a combined analysis to confirm that the studied taxa are genetically arranged following the geographic speciation centers. After verifying this premise, we next examined whether the current diversity of species of Phalolepis in the four Mediterranean peninsulas agreed with a MGH model. The first condition (the presence of lowland, montane and alpine zones) is met in the four speciation centers, and thereafter, we focused on the second and third conditions of this model. To explore the second condition, i.e., the presence of climatic fluctuations that enabled mountains to act as "species pump", we used the recently published Climatic Stability Index (CSI) [22] to examine whether there is a correlation between climatic instability and the impetus of allopatric speciation. Finally, we examined the components of the third condition, namely high relief terrain (ruggedness) and environmental gradients (climate niche breadth), which were also tested against the rate of allopatric speciation by region.

Plant Material
The total number of individuals studied was 1487 from 56 populations belonging to 21 species, all of them diploid ( Figure 1). The number of species by region was as follows: Iberian Peninsula, two species (C. alba L. and C. costae Willk.); Italian Peninsula + Northern Balkans, six species (C. aspromontana Brullo, Scelsi & Spamp., C. ionica Brullo, C. nobilis (Groves) Brullo . Centaurea deusta and C. alba are species with a wide distribution compared with the rest, which are narrow endemics. The voucher information is shown in Supplemental  Table S1. The Italian Peninsula and the Northern Balkans were merged in a single unit for this analysis because of the evidence of recent connections between both regions [23].

Genetic Structure and Allopatric Speciation Parameters
For the analysis of the genetic structure, the clustering software STRUCTURE v. 2.3.4 [26], which is based on a Bayesian approach, was used. In total, 1487 individuals (from 56 populations from the four regions) were analyzed. After a preliminary run with the number of clusters (K) set from 1 to 56 and 10 iterations per K, we restricted K to 1 to 37 and increased the number iterations per K to 20, assuming an admixture model with correlated allele frequencies and no "locprior". The burn-in length and the number of Markov chain Monte Carlo (MCMC) replications were set to 10 5 and 10 6 , respectively. The most likely value of K was determined in two ways: by choosing the smallest K according to the log-probability of the data (ln Pr(X|K)) where the values reached a plateau [27], and by the ∆K statistical approach of [28] implemented in the software STRUCTURE HARVESTER v. 0.6.94 [29]. Finally, the programs Clumpp v. 1.1.2 [30] and Distruct v. 1.1 [31] were used in order to obtain a clear and interpretable image of this analysis. Additionally, the software GenAlEx v. 6.5 [32] was used to perform a principal coordinate analysis (PCoA) at the population level, based on the codominant genotype distances.

Genetic Structure and Allopatric Speciation Parameters
For the analysis of the genetic structure, the clustering software STRUCTURE v. 2.3.4 [26], which is based on a Bayesian approach, was used. In total, 1487 individuals (from 56 populations from the four regions) were analyzed. After a preliminary run with the number of clusters (K) set from 1 to 56 and 10 iterations per K, we restricted K to 1 to 37 and increased the number iterations per K to 20, assuming an admixture model with correlated allele frequencies and no "locprior". The burn-in length and the number of Markov chain Monte Carlo (MCMC) replications were set to 10 5 and 10 6 , respectively. The most likely value of K was determined in two ways: by choosing the smallest K according to the log-probability of the data (ln Pr(X|K)) where the values reached a plateau [27], and by the ΔK statistical approach of [28] implemented in the software STRUCTURE HARVESTER v. 0.6.94 [29]. Finally, the programs Clumpp v. 1.1.2 [30] and Distruct v. 1.1 [31] were used in order to obtain a clear and interpretable image of this analysis. Additionally, the software GenAlEx v. 6.5 [32] was used to perform a principal coordinate analysis (PCoA) at the population level, based on the codominant genotype distances.
To infer the degree of allopatric speciation in Centaurea subsect. Phalolepis, we compared the levels of genetic diversity, genetic diversification, recent gene flow and past gene flow for the four speciation areas [10][11][12][13]. Specifically, we calculated the expected heterozygosity (He) as the best estimator of genetic diversity [33], using GenAlEx. For genetic differentiation between populations, we chose FST, which was calculated using FreeNA [34] and/or Arlequin v. 3.11 [35]. For recent gene flow, we used recent (last onethree generations) migration rates (m) estimated using the program BayesAss v. 1.3 and To infer the degree of allopatric speciation in Centaurea subsect. Phalolepis, we compared the levels of genetic diversity, genetic diversification, recent gene flow and past gene flow for the four speciation areas [10][11][12][13]. Specifically, we calculated the expected heterozygosity (H e ) as the best estimator of genetic diversity [33], using GenAlEx. For genetic differentiation between populations, we chose F ST, which was calculated using FreeNA [34] and/or Arlequin v. 3.11 [35]. For recent gene flow, we used recent (last onethree generations) migration rates (m) estimated using the program BayesAss v. 1.3 and v. 3.0 [36]. For older rates of gene flow, we calculated historical mutation-scaled migration rates (M), originally obtained with the software MIGRATE-N [37]. The effective number of migrants per generation (Nm) was estimated using the formula 4Nm = ΘM [38]. For more explanations and descriptions of the genetic parameters and the methods used to estimate them, please refer to the four cited studies.

Testing the MGH Model
To check whether the second and third conditions (climatic fluctuations and highrelief terrain with large environmental gradients, respectively) of the MGH model were responsible for the higher degree of allopatric speciation of the eastern compared with western speciation centers, we calculated the values for adequate surrogates of three conditions and we searched for statistical differences among the four geographic regions. For the second condition, we chose the Climatic Stability Index (CSI) [22]. The "CSIpast" uses 14 of the 19 bioclimatic variables of the WorldClim database (https://www. worldclim.org, accessed 20 November 2020) for the 12 periods of PaleoClim [39] that represent warm/cold cycles from the Pliocene (ca. 3.3 Ma) to the present; it is available at a 2.5 arc-min (~5 km) resolution covering the whole planet, with values between 0 (completely stable areas) and 1 (the most unstable ones). We believe that the CSI is particularly suited Plants 2023, 12, 11 5 of 15 for our species system because speciation events within Centaurea subsection Phalolepis took place from the mid-Pliocene to the present [40,41] within the timespan covered by the CSI. The CSI was clipped to an area delimited by 50 • N, 12 • W, 33 • N and 42 • E in R using the package Raster v.3.5-15 [42].
For the third condition, we chose one variable for each of the two subconditions. For representing high-relief terrain, we chose the terrain ruggedness index (TRI), which is the mean of the absolute differences in elevation between a focal cell and its eight surrounding cells; flat areas have a value of zero. The TRI was obtained from [43], and it is available at several resolutions (we used a 2.5 arc-min resolution). For the second subcondition (environmental gradients), we generated an ad hoc variable, the Climate Niche Breadth (CNB), which was estimated by using all 19 bioclimatic variables (at a 30 arc-sec resolution) of the WorldClim database (accessed on 28 February 2022). The CNB was calculated by using the R package base v.4.1.1 [44] as follows. First, we calculated the standard deviation (SD) for each bioclimatic variable selected (see below) from the 25 pixels at 30 arc-sec resolution that were included within a 2.5 arc-min pixel. We then normalized all the values for each 2.5 arc-min pixel between 0 and 1. Finally, we summed all the SD values for each selected variable (see below). Higher CNB values represented broad (climate) niches, whereas lower values indicated narrow niches. Both the TRI and CNB were clipped to the same study area as the CSI.
As there could be problems associated with the collinearity of bioclimatic variables, we performed a variable selection process to reduce the autocorrelation effect for estimating both the CSI and CNB. With the clipped variables confined to the study area, a Pearson correlation analysis was performed to identify the possibly highly correlated variable sets using the R package Stats v.4.1.1 [44]. Next, a principal component analysis (PCA) was conducted by using the R package FactoMineR v.2.4 [45] to reduce the dimensions of the explanatory variables. The contribution of the variables to each of the PCA axes was obtained with the R package factoextra v.1.0.7. [46]. We retained the first three principal components (PCs) that accounted for >80% of the variance in the original data [47]. The three variables that contributed the most to explaining variance in the first three PCs were conserved; of these, we discarded those with the least percentage of contribution to that axis from highly correlated sets (Pearson's r >|0.9|). For the CSI, six final variables were used for further analysis: bio8 (the mean temperature of the wettest quarter), bio11 (the mean temperature of the coldest quarter), bio16 (the precipitation of the wettest quarter), bio17 (the precipitation of the driest quarter), bio18 (the precipitation of the warmest quarter) and bio19 (the precipitation of the coldest quarter). For the CNB, the variables selected, in addition to bio8, bio17 and bio19, were bio5 (the maximum temperature of the warmest month), bio6 (the minimum temperature of the coldest month) and bio15 (seasonality of precipitation).
The values for the three variables (CSI, TRI and CNB) for each of the speciation centers were calculated by three different strategies: (1) using only the species occurrence points, i.e., the values for the three variables were extracted for each occurrence with the "extract" function of the raster package in R and (2) delimiting 10 km buffer areas around each occurrence. Comparisons among the speciation centers were statistically tested. First, the Shapiro-Wilk test was performed to check the normality of the data; as they were not normally distributed (the p value was less than an alpha level of 0.05), non-parametric Wilcoxon-Mann-Whitney analysis was conducted. The analyses were performed using the shapiro.test and wilcox.test functions of the package Stats v. 4.1.1 in R.

Genetic Structure and Allopatric Speciation Parameters
According to the approach of [28] for calculating the best number of clusters (K), the most likely value of K was 4 (Supplementary Figure S1), in agreement with the four biogeographic areas of subsection Phalolepis: Greece (cream color in Figure 2), Turkey (blue), the Iberian Peninsula (red) and Italy and the Northern Balkans (green), although this last region showed a high degree of admixture, specifically in the Northern Balkans ( Figure 2). In the principal coordinate analysis (PCoA) at population level (Figure 3), the resolution was good as far as the Greek and Turkish populations were concerned. However, the Balkan and Italian populations ION2, DEU3, DEU4, DEU5, DEU6 and DEU7 showed a clear overlapping zone with the Iberian populations in Figure 3; in these six populations, the presence of the Iberian genetic cluster (indicated by the red color in Figure 2) was evident in the STRUCTURE analysis.

Genetic Structure and Allopatric Speciation Parameters
According to the approach of [28] for calculating the best number of clusters (K), the most likely value of K was 4 (Supplementary Figure S1), in agreement with the four biogeographic areas of subsection Phalolepis: Greece (cream color in Figure 2), Turkey (blue), the Iberian Peninsula (red) and Italy and the Northern Balkans (green), although this last region showed a high degree of admixture, specifically in the Northern Balkans ( Figure 2). In the principal coordinate analysis (PCoA) at population level (Figure 3), the resolution was good as far as the Greek and Turkish populations were concerned. However, the Balkan and Italian populations ION2, DEU3, DEU4, DEU5, DEU6 and DEU7 showed a clear overlapping zone with the Iberian populations in Figure 3; in these six populations, the presence of the Iberian genetic cluster (indicated by the red color in Figure 2) was evident in the STRUCTURE analysis.   Regarding the genetic diversity, the Iberian Peninsula and Northern Balkans-Italy were the areas with higher genetic variation (measured as the expected heterozygosity He) and both past and recent gene flow (Table 1). In particular, the Iberian Peninsula had clearly much higher values (about one order of magnitude higher) for both recent and past gene flow than the other areas. For genetic differentiation (FST), the values from our previous analyses were the highest in Greece and lowest in the Iberian Peninsula (FST = 0.243 and 0.176, respectively; see Table 1).  Regarding the genetic diversity, the Iberian Peninsula and Northern Balkans-Italy were the areas with higher genetic variation (measured as the expected heterozygosity H e ) and both past and recent gene flow (Table 1). In particular, the Iberian Peninsula had clearly much higher values (about one order of magnitude higher) for both recent and past gene flow than the other areas. For genetic differentiation (F ST ), the values from our previous analyses were the highest in Greece and lowest in the Iberian Peninsula (F ST = 0.243 and 0.176, respectively; see Table 1).

Testing the Mountain-Geobiodiversity Hypothesis (MGH) Model
We checked the impact of climatic fluctuations and high-relief terrain with large environmental gradients as the second and third conditions of the MGH using the Climate Stability Index (CSI), the Terrain Ruggedness Index (TRI) and the Climate Niche Breadth (CNB). For both approaches (species occurrences and 10 km buffer areas), we found statistical differences for all comparisons regarding the Climatic Stability Index (CSI), with the only exception of Northern Balkans-Italy vs. Greece for the dataset of occurrences (Figure 4 and Supplementary Table S2). The largest values of CSI corresponded to Northern Balkans-Italy, followed by the Iberian Peninsula, Greece and Turkey. For the TRI, all comparisons were statistically significant for both occurrences and 10 km buffers; Greece was the most rugged area, followed by Turkey, Northern Balkans-Italy and the Iberian Peninsula ( Figure 4 and Supplementary Table S2). Regarding the CNB, Greece showed the widest niche breadth, followed by Turkey, Northern Balkans-Italy and the Iberian Peninsula; however, while we detected statistical differences for all comparison pairs when taking 10 km buffers into account, for the occurrences, differences were only found when we compared the Iberian Peninsula vs. Northern Balkans-Italy and Greece (Figure 4 and Supplementary Table S2).  Table S2). The largest values of CSI corresponded to Northern Balkans-Italy, followed by the Iberian Peninsula, Greece and Turkey. For the TRI, all comparisons were statistically significant for both occurrences and 10 km buffers; Greece was the most rugged area, followed by Turkey, Northern Balkans-Italy and the Iberian Peninsula (Figure 4 and Supplementary Table S2). Regarding the CNB, Greece showed the widest niche breadth, followed by Turkey, Northern Balkans-Italy and the Iberian Peninsula; however, while we detected statistical differences for all comparison pairs when taking 10 km buffers into account, for the occurrences, differences were only found when we compared the Iberian Peninsula vs. Northern Balkans-Italy and Greece ( Figure 4 and Supplementary Table S2).

Genetic Parameters and Speciation
The genetic structure resulting from our analyses reflects the presence of four geographic groups. However, there are some discordant populations in the Northern Balkans-Italy group, which show high levels of admixture with a predominance of the Iberian red genotype (Figure 2). The same Northern Balkans-Italy populations

Genetic Parameters and Speciation
The genetic structure resulting from our analyses reflects the presence of four geographic groups. However, there are some discordant populations in the Northern Balkans-Italy group, which show high levels of admixture with a predominance of the Iberian red genotype (Figure 2). The same Northern Balkans-Italy populations overlapped with some Iberian populations in the PCoA (Figure 3). The similarity of both genetic groups could be the result of hybridization of these populations with other species from subsect. Centaurea. In support of this hypothesis, we detected very high levels of admixture in the Northern Balkan-Italian populations [12], and hybridization was verified in the Iberian C. costae [13]. Genetic diversity is higher in the Iberian and Northern Balkans-Italy populations compared with the populations of the other speciation centers ( Table 1). The high levels of genetic diversity in Centaurea could also be the result of hybridization followed by the repeated backcrossing of interspecific hybrids with one of the parents, as suggested by [48] for populations of Centaurea podospermifolia Loscos & J. Pardo. Certainly, increases in genetic diversity have frequently been reported in hybridizing populations in plants [49]. Both historical and recent interspecific gene flow in the Iberian Peninsula are extremely high, up to 10 times higher than in the other regions (Table 1), which support the extent of introgression within Iberian populations. High levels of gene flow are important to maintain high levels of genetic diversity within species, as alleles are less likely to be eliminated; due to high levels of gene flow, alleles tend to be shared in several populations, thus minimizing the effects of genetic drift. Only in the Iberian Peninsula, the historical levels of gene flow are theoretically capable of counteracting the effects of genetic drift (Nm > 1; [50]).
The Turkish and Greek populations, on the contrary, show the lowest levels of both historical and current gene flow, the lowest values of H e and the highest values of F ST , with the exception of Turkey for F ST , thus confirming that the species in these two speciation centers are strongly isolated (Table 1). According to the genetic data, we can conclude that speciation on the Phalolepis group follows two different models: allopatry by isolation in rugged mountains, exemplified by the cases of Greece and Turkey, and hybrid speciation as evidenced for the Iberian Peninsula. Speciation in the Northern Balkans-Italy region is closer to the Iberian model, even though gene flow is lower (Table 1).

Mountains and Speciation: Is the MGH Model Applicable to Our Case Study?
Although the elegant MGH model was formulated for explaining the factors accounting for the high biodiversity richness in the Tibeto-Himalayan region (THR), it was tested shortly after on 16 mountain systems around the world harboring important plant radiations [4]. With the help of general linear models (GLMs), the authors concluded that their model was supported, as the combination of three parameters was statistically significant for explaining plant species richness in mountain ranges: (1) net primary productivity (NPP), (2) elevation range, and (3) the change in the mean annual temperature (∆T) between the LGM (Last Glacial Maximum, ca. 20,000 years BP) and today (the latter was taken as a rough proxy for climatic fluctuations for the whole Quaternary Period, i.e., the last 2.6 My). In addition, the authors also found that the range of elevation was positively correlated with a geodiversity index. However, ∆T showed a negative correlation with species richness; the authors [4] concluded that "mountain systems characterized by small elevation ranges and strong modifications of temperature [ . . . ] appear to harbor less radiations", and most radiations were identified in high-elevation mountains and weaker temperature oscillations. These patterns do not deny the second condition of the MGH hypothesis (the "species pump" effect associated with climatic instability) because, as noted by [4], most shifts in diversification rates of the radiations examined by these authors occurred much later than mountain-building processes. For the examples of the Tibetan-Hengduan Region and the Andes, this means that many radiations took place a few My after the orogenic movements of the Miocene, coinciding with the climatic changes of the Pliocene and Pleistocene Epochs. Alternatively, the negative correlation of species richness with ∆T could be due to the inability of this parameter (which measures the temperature differences between ca. 21,000 years ago and the present) to infer the temperature variations through the Pleistocene and Pliocene Epochs, when the vast majority of the 82 radiations examined by [4] took place; the cyclic temperature oscillations (∆T) progressively increased from just about 1 • C in the early Pliocene to >6 • C in the late Pleistocene [51,52]. The recently developed CSI [22] uses 12 periods from the mid-Pliocene (3.3 My) to the present to estimate climate stability/instability. Thus, it could serve as a better option than ∆T LGM-present for assessing the effect of climatic instability as the trigger of the species pump.
Our results for Centaurea in the Mediterranean mountains, which suggest a positive correlation between mountain ruggedness and speciation, fit well with the effects of ∆T LGM-present observed by [4]. Although the CSI index used here combines six climatic variables instead of a single one, the regions with the maximum species diversity of subsect. Phalolepis (Turkey and Greece) have been the most climatically stable since 3.3 Mya, whereas the most climatically unstable regions (Iberian Peninsula, Italy-Northern Balkans) are the poorest in terms of species richness. The third condition of the MGH hypothesis ("high-relief terrain with environmental gradients") is clearly met in the mountain systems of Turkey and Greece: these two regions showed the highest values for both TRI and CNB.
The topographical and climatic characteristics of the Iberian Peninsula and, at a smaller scale, Northern Balkans-Italy made them ideal for hybridogenic speciation; the combination of the strong climatic oscillations and the gentle topography forced and allowed, respectively, latitudinal shifts in species (Scenario 2 in Figure 5). Large migration movements caused species contact and favored hybridization; the special topography of the central part of the Iberian Peninsula would have enabled migrations of hundreds of km, which would have ensured the integrity of, e.g., Centaurea alba [13]. The very low number of species of subsect. Phalolepis in this region is not unexpected and is likely due to the limited speciation processes via hybridization coupled with putative extinction processes; some species, particularly narrow endemics, would have not had the chances for altitudinal migrations in the low-elevation Meseta ( Figure 5). movements caused species contact and favored hybridization; the special topography of the central part of the Iberian Peninsula would have enabled migrations of hundreds of km, which would have ensured the integrity of, e.g., Centaurea alba [13]. The very low number of species of subsect. Phalolepis in this region is not unexpected and is likely due to the limited speciation processes via hybridization coupled with putative extinction processes; some species, particularly narrow endemics, would have not had the chances for altitudinal migrations in the low-elevation Meseta ( Figure 5).  Greece and Turkey, in contrast, became relatively important speciation centers of subsect. Phalolepis, with events of allopatric speciation shown by the genetic markers, likely the consequence of a very rugged terrain providing geographical barriers. The speciation of Centaurea does not constitute a radiation in any way because the nuclei of maximum diversity in Anatolia and Greece comprise no more than 30 species, the mountains have only moderate altitudes of generally 2000-3000 m, and the climate was probably not unstable enough. Nevertheless, within the context of the Mediterranean Basin, both the Greek and Turkish mountain ranges may superficially resemble the THR and the Andes [4]. These two very high cordilleras are examples of "mountain systems with the largest elevational ranges and stronger overlap between today's and LGM temperature profiles and are also those where most radiations were identified" [4]. Rugged mountains with large altitudinal gradients allow species to survive climate changes by moving up and down; if the climate oscillations are not large, the required altitudinal shifts are not so extensive, which would help populations to keep populations large enough to avoid local extinction (Scenario 3 in Figure 5). If high ecological gradients are added to the equation (herein measured though the CNB), then plant species could easily find favorable pockets for survival (e.g., "microrefugia" sensu [53]) while stimulating ecological speciation [54,55]. Allopatric and ecological speciation are not mutually exclusive; thus, they can contribute to the patterns of high species richness in topographically and ecologically complex mountains.
The patterns of Centaurea subsect. Phalolepis in the Turkish and Greek mountain ranges greatly resemble those of Scenario 3 in Figure 5: the dual role of mountains as species cradles due to allopatric speciation, and as species museums due to the occurrence of (micro)refugia. This feature of many mountains has undoubtedly attracted much interest on a conservation basis [56][57][58], as their efficient preservation would allow for the protection of those ecological and evolutionary processes that create and maintain biodiversity due to the inherent topographic and environmental heterogeneity.

Climate Niche Breadth (CNB): A New Index for Representing Climate Niche Breadth at the Global Scale
Here, we have implemented a new parameter that has allowed us to test environmental gradients against the rate of allopatric speciation in Centaurea subsect. Phalolepis. Although the CNB only includes climatic aspects, it has the paramount advantage that data are available for the whole Earth (the bioclimatic variables of WorldClim); therefore, users could select a desired subset of the climatic variables (or all 19 combined ones) and create their "own" CNB for a given region (or the whole planet). Through Figshare (https: //doi.org/10.6084/m9.figshare.20863528), researchers can access two sets of data records: (1) SD-based maps of individual bioclimatic variables, which contain a set of raster layers (19, one for each variable); (2) a SD-based map for the 19 bioclimatic variables combined; and (3) the scripts that allow them to combine the raster layers according to the user's preferences to generate a customized CNB.
As expected, there is a large positive correlation between TRI and CNB: r = 0.796 for our study area. By definition, rugged terrains provide a wide array of climatic conditions; thus, it is not rare that rugged mountains have traditionally been linked to glacial refugia, as species were able to find local microclimates and thus persisted [59,60]. Despite this relationship between TRI and CNB, the r values are well below 1, which indicates that the ruggedness alone cannot fully explain the breadth of climatic niches. This correlation is reduced on a global level (r = 0.634); if we compare the global maps of TRI and CNB ( Figure 6), there are some important differences. For instance, the Colombian Andes, which is one the regions of the world with the highest levels of CNB, has only moderate values for TRI, and the same is applicable to the Alps. Conversely, some areas with high/moderate values of TRI show moderate/low values for CNB, particularly in central, east (including the Russian Far East) and southeast Asia ( Figure 6). The differences can be very high at a more local scale; examples include some parts of northern Queensland (Australia) and Pantepui in the Guiana Region, with the values of CNB much higher than those of TRI ( Figure 6). In the case of the Pantepui region, this discrepancy is easy to explain: the topography is relatively simple, and the landscape is characterized by the tabular mountains or "tepuis", which are typically 2000-2600 m high, emerging from the lowland savannas, which are generally 100-400 m high [61], but these large vertical gradients offer a wide array of microclimates. (Australia) and Pantepui in the Guiana Region, with the values of CNB much higher than those of TRI ( Figure 6). In the case of the Pantepui region, this discrepancy is easy to explain: the topography is relatively simple, and the landscape is characterized by the tabular mountains or "tepuis", which are typically 2000-2600 m high, emerging from the lowland savannas, which are generally 100-400 m high [61], but these large vertical gradients offer a wide array of microclimates.

Concluding Remarks
Our analyses of Phalolepis highlighted the crucial role of the mountains in both the genesis and the preservation of species diversity. Ruggedness emerged as the dominant factor, but all the conditions indicated by [3] for the Mountain-Geobiodiveristy Hypothesis, namely the presence of lowland, montane and alpine zones; climatic oscillations that enable mountains to act as species pumps; and high-relief terrain with environmental gradients, must be present for true radiation. In the Iberian Peninsula, with the lowest ruggedness and the presence of a high plateau connecting mountain ranges, the diversity was very low and speciation was hybridogenic. In the other Mediterranean mountains, high ruggedness favored allopatric speciation; however, climatic stability was high, which would have impeded true radiation and, at the same time, reinforced their role as refugia during periods of glaciation. This pattern of species evolution stresses the need for special protection of the montane ranges in the Mediterranean Region.

Concluding Remarks
Our analyses of Phalolepis highlighted the crucial role of the mountains in both the genesis and the preservation of species diversity. Ruggedness emerged as the dominant factor, but all the conditions indicated by [3] for the Mountain-Geobiodiveristy Hypothesis, namely the presence of lowland, montane and alpine zones; climatic oscillations that enable mountains to act as species pumps; and high-relief terrain with environmental gradients, must be present for true radiation. In the Iberian Peninsula, with the lowest ruggedness and the presence of a high plateau connecting mountain ranges, the diversity was very low and speciation was hybridogenic. In the other Mediterranean mountains, high ruggedness favored allopatric speciation; however, climatic stability was high, which would have impeded true radiation and, at the same time, reinforced their role as refugia during periods of glaciation. This pattern of species evolution stresses the need for special protection of the montane ranges in the Mediterranean Region.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/plants12010011/s1, Supplementary Table S1: Voucher information. Supplementary Figure S1: Results of the Evanno test [28] for the best K. Supplementary  Table S2: p values for the comparisons among the four speciation centers regarding the Climate Stability Index (CSI), the Terrain Ruggedness Index (TRI) and the Climate Niche Breadth (CNB).  Data Availability Statement: All the data used in this paper are already available; please refer to [10][11][12][13].

Conflicts of Interest:
The authors declare no conflict of interest.